Stokes I Simultaneous Image and Instrument Modeling

In this tutorial, we will create a preliminary reconstruction of the 2017 M87 data on April 6 by simultaneously creating an image and model for the instrument. By instrument model, we mean something akin to self-calibration in traditional VLBI imaging terminology. However, unlike traditional self-cal, we will at each point in our parameter space effectively explore the possible self-cal solutions. This will allow us to constrain and marginalize over the instrument effects, such as time variable gains.

To get started we load Comrade.

using Comrade
  Activating project at `~/work/Comrade.jl/Comrade.jl/examples`
using Pyehtim
using LinearAlgebra

For reproducibility we use a stable random number genreator

using StableRNGs
rng = StableRNG(123)
StableRNGs.LehmerRNG(state=0x000000000000000000000000000000f7)

Load the Data

To download the data visit https://doi.org/10.25739/g85n-f134 First we will load our data:

obs = ehtim.obsdata.load_uvfits(joinpath(dirname(pathof(Comrade)), "..", "examples", "SR1_M87_2017_096_hi_hops_netcal_StokesI.uvfits"))
Python: <ehtim.obsdata.Obsdata object at 0x7f5ce0f29f90>

Now we do some minor preprocessing:

  • Scan average the data since the data have been preprocessed so that the gain phases coherent.
  • Add 1% systematic noise to deal with calibration issues that cause 1% non-closing errors.
obs = scan_average(obs.add_fractional_noise(0.02))
Python: <ehtim.obsdata.Obsdata object at 0x7f5cddfd4e50>

Now we extract our complex visibilities.

dvis = extract_table(obs, ComplexVisibilities())
EHTObservation{Float64,Comrade.EHTVisibilityDatum{Float64}, ...}
  source: M87
  mjd: 57849
  frequency: 2.29070703125e11
  bandwidth: 1.856e9
  stations: [:AA, :AP, :AZ, :JC, :LM, :PV, :SM]
  nsamples: 284

##Building the Model/Posterior

Now, we must build our intensity/visibility model. That is, the model that takes in a named tuple of parameters and perhaps some metadata required to construct the model. For our model, we will use a raster or ContinuousImage for our image model. The model is given below:

function sky(θ, metadata)
    (;fg, c, σimg) = θ
    (;ftot, K, meanpr, grid, cache) = metadata
    # Transform to the log-ratio pixel fluxes
    cp = meanpr .+ σimg.*c.params
    # Transform to image space
    rast = (ftot*(1-fg))*K(to_simplex(CenteredLR(), cp))
    img = IntensityMap(rast, grid)
    m = ContinuousImage(img, cache)
    # Add a large-scale gaussian to deal with the over-resolved mas flux
    g = modify(Gaussian(), Stretch(μas2rad(250.0), μas2rad(250.0)), Renormalize(ftot*fg))
    return m + g
end
sky (generic function with 1 method)

Unlike other imaging examples (e.g., Imaging a Black Hole using only Closure Quantities) we also need to include a model for the instrument, i.e., gains as well. The gains will be broken into two components

  • Gain amplitudes which are typically known to 10-20%, except for LMT, which has amplitudes closer to 50-100%.
  • Gain phases which are more difficult to constrain and can shift rapidly.
function instrument(θ, metadata)
    (; lgamp, gphase) = θ
    (; gcache, gcachep) = metadata
    # Now form our instrument model
    gvis = exp.(lgamp)
    gphase = exp.(1im.*gphase)
    jgamp = jonesStokes(gvis, gcache)
    jgphase = jonesStokes(gphase, gcachep)
    return JonesModel(jgamp*jgphase)
end
instrument (generic function with 1 method)

The model construction is very similar to Imaging a Black Hole using only Closure Quantities, except we include a large scale gaussian since we want to model the zero baselines. For more information about the image model please read the closure-only example. Let's discuss the instrument model Comrade.JonesModel. Thanks to the EHT pre-calibration, the gains are stable over scans. Therefore, we can model the gains on a scan-by-scan basis. To form the instrument model, we need our

  1. Our (log) gain amplitudes and phases are given below by lgamp and gphase
  2. Our function or cache that maps the gains from a list to the stations they impact gcache.
  3. The set of Comrade.JonesPairs produced by jonesStokes

These three ingredients then specify our instrument model. The instrument model can then be combined with our image model cimg to form the total JonesModel.

Now, let's set up our image model. The EHT's nominal resolution is 20-25 μas. Additionally, the EHT is not very sensitive to a larger field of view. Typically 60-80 μas is enough to describe the compact flux of M87. Given this, we only need to use a small number of pixels to describe our image.

npix = 32
fovx = μas2rad(150.0)
fovy = μas2rad(150.0)
7.27220521664304e-10

Now let's form our cache's. First, we have our usual image cache which is needed to numerically compute the visibilities.

grid = imagepixels(fovx, fovy, npix, npix)
buffer = IntensityMap(zeros(npix, npix), grid)
cache = create_cache(NFFTAlg(dvis), buffer, BSplinePulse{3}())
VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}(VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}(NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}(1, 4, 2.0, :kaiser_bessel, AbstractNFFTs.POLYNOMIAL, true, false, true, 0x00000000), [1.2065151329523807e9 819926.8691406249 … 4.472903950222221e9 -15688.352240668402; -3.8185738483809524e9 -1.7123166484375e6 … -6.2063795222222224e7 119336.55902777778]), NFFTPlan with 284 sampling points for an input array of size(32, 32) and an output array of size(284,) with dims 1:2, ComplexF64[0.9304270810105981 - 0.1755525366839422im, 0.9999999857224657 - 6.37119066797087e-5im, 0.6968566352375541 - 0.5221791050707502im, 0.9304860627318433 + 0.17550202040945986im, 0.8028821861479044 + 0.3950138433443315im, 0.696909133455242 + 0.522149977275377im, 0.9999999859666728 - 6.454986502401626e-5im, 0.690100908057055 - 0.5259755617707915im, 0.9317353097306306 - 0.1646343527093735im, 0.9317937328057632 + 0.1645830409592267im  …  0.8346090489457091 + 0.1795789988845693im, 0.834610800886739 + 0.17958571977595708im, 0.9181402138641586 + 0.09206131734480663im, 0.9120864418553256 + 0.2942607208902332im, 0.9120838428437319 + 0.29426742718619836im, 0.9559925804493471 + 0.09847981391058878im, 0.8881131676530709 + 0.289296040490046im, 0.9924358388500316 + 0.0028007013218525043im, 0.8881114735852544 + 0.28930284794286826im, 0.9999999999233905 + 7.3999334332420275e-6im], BSplinePulse{3}(), [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])

Second, we now construct our instrument model cache. This tells us how to map from the gains to the model visibilities. However, to construct this map, we also need to specify the observation segmentation over which we expect the gains to change. This is specified in the second argument to jonescache, and currently, there are two options

  • FixedSeg(val): Fixes the corruption to the value val for all time. This is usefule for reference stations
  • ScanSeg(): which forces the corruptions to only change from scan-to-scan
  • TrackSeg(): which forces the corruptions to be constant over a night's observation

For this work, we use the scan segmentation for the gain amplitudes since that is roughly the timescale we expect them to vary. For the phases we use a station specific scheme where we set AA to be fixed to unit gain because it will function as a reference station.

gcache = jonescache(dvis, ScanSeg())
gcachep = jonescache(dvis, ScanSeg(); autoref=SEFDReference((complex(1.0))))

using VLBIImagePriors

Now we need to specify our image prior. For this work we will use a Gaussian Markov Random field prior Since we are using a Gaussian Markov random field prior we need to first specify our mean image. This behaves somewhat similary to a entropy regularizer in that it will start with an initial guess for the image structure. For this tutorial we will use a a symmetric Gaussian with a FWHM of 50 μas

fwhmfac = 2*sqrt(2*log(2))
mpr = modify(Gaussian(), Stretch(μas2rad(50.0)./fwhmfac))
imgpr = intensitymap(mpr, grid)
2-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓   X ∈ 32-element LinRange{Float64,...}
→   Y ∈ 32-element LinRange{Float64,...}
And data, 32×32 NamedDimsArray(::Matrix{Float64}, (:X, :Y)):
                 (-3.52247e-10)  …  (3.29522e-10)  (3.52247e-10)
 (-3.52247e-10)     6.37537e-8        1.32434e-7     6.37537e-8
 (-3.29522e-10)     1.32434e-7        2.751e-7       1.32434e-7
 (-3.06796e-10)     2.62014e-7        5.44273e-7     2.62014e-7
 (-2.84071e-10)     4.93724e-7        1.0256e-6      4.93724e-7
 (-2.61345e-10)     8.86092e-7   …    1.84065e-6     8.86092e-7
 (-2.38619e-10)     1.51463e-6        3.14629e-6     1.51463e-6
    ⋮                            ⋱    ⋮            
  (2.15894e-10)     2.46586e-6        5.12225e-6     2.46586e-6
  (2.38619e-10)     1.51463e-6        3.14629e-6     1.51463e-6
  (2.61345e-10)     8.86092e-7   …    1.84065e-6     8.86092e-7
  (2.84071e-10)     4.93724e-7        1.0256e-6      4.93724e-7
  (3.06796e-10)     2.62014e-7        5.44273e-7     2.62014e-7
  (3.29522e-10)     1.32434e-7        2.751e-7       1.32434e-7
  (3.52247e-10)     6.37537e-8        1.32434e-7     6.37537e-8

Now since we are actually modeling our image on the simplex we need to ensure that our mean image has unit flux

imgpr ./= flux(imgpr)
2-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓   X ∈ 32-element LinRange{Float64,...}
→   Y ∈ 32-element LinRange{Float64,...}
And data, 32×32 Matrix{Float64}:
                 (-3.52247e-10)  …  (3.29522e-10)  (3.52247e-10)
 (-3.52247e-10)     6.38049e-8        1.3254e-7      6.38049e-8
 (-3.29522e-10)     1.3254e-7         2.75321e-7     1.3254e-7
 (-3.06796e-10)     2.62224e-7        5.4471e-7      2.62224e-7
 (-2.84071e-10)     4.94121e-7        1.02642e-6     4.94121e-7
 (-2.61345e-10)     8.86803e-7   …    1.84213e-6     8.86803e-7
 (-2.38619e-10)     1.51585e-6        3.14882e-6     1.51585e-6
    ⋮                            ⋱    ⋮            
  (2.15894e-10)     2.46784e-6        5.12636e-6     2.46784e-6
  (2.38619e-10)     1.51585e-6        3.14882e-6     1.51585e-6
  (2.61345e-10)     8.86803e-7   …    1.84213e-6     8.86803e-7
  (2.84071e-10)     4.94121e-7        1.02642e-6     4.94121e-7
  (3.06796e-10)     2.62224e-7        5.4471e-7      2.62224e-7
  (3.29522e-10)     1.3254e-7         2.75321e-7     1.3254e-7
  (3.52247e-10)     6.38049e-8        1.3254e-7      6.38049e-8

and since our prior is not on the simplex we need to convert it to unconstrained or real space.

meanpr = to_real(CenteredLR(), Comrade.baseimage(imgpr))
32×32 Matrix{Float64}:
 -7.55422  -6.82317  -6.14085  -5.50727   …  -6.14085  -6.82317  -7.55422
 -6.82317  -6.09211  -5.4098   -4.77622      -5.4098   -6.09211  -6.82317
 -6.14085  -5.4098   -4.72748  -4.0939       -4.72748  -5.4098   -6.14085
 -5.50727  -4.77622  -4.0939   -3.46032      -4.0939   -4.77622  -5.50727
 -4.92243  -4.19137  -3.50906  -2.87548      -3.50906  -4.19137  -4.92243
 -4.38632  -3.65527  -2.97295  -2.33937   …  -2.97295  -3.65527  -4.38632
 -3.89895  -3.1679   -2.48558  -1.852        -2.48558  -3.1679   -3.89895
 -3.46032  -2.72927  -2.04695  -1.41337      -2.04695  -2.72927  -3.46032
 -3.07043  -2.33937  -1.65705  -1.02348      -1.65705  -2.33937  -3.07043
 -2.72927  -1.99821  -1.3159   -0.682317     -1.3159   -1.99821  -2.72927
  ⋮                                       ⋱             ⋮        
 -3.07043  -2.33937  -1.65705  -1.02348      -1.65705  -2.33937  -3.07043
 -3.46032  -2.72927  -2.04695  -1.41337      -2.04695  -2.72927  -3.46032
 -3.89895  -3.1679   -2.48558  -1.852     …  -2.48558  -3.1679   -3.89895
 -4.38632  -3.65527  -2.97295  -2.33937      -2.97295  -3.65527  -4.38632
 -4.92243  -4.19137  -3.50906  -2.87548      -3.50906  -4.19137  -4.92243
 -5.50727  -4.77622  -4.0939   -3.46032      -4.0939   -4.77622  -5.50727
 -6.14085  -5.4098   -4.72748  -4.0939       -4.72748  -5.4098   -6.14085
 -6.82317  -6.09211  -5.4098   -4.77622   …  -5.4098   -6.09211  -6.82317
 -7.55422  -6.82317  -6.14085  -5.50727      -6.14085  -6.82317  -7.55422

Now we can form our metadata we need to fully define our model.

metadata = (;ftot=1.1, K=CenterImage(imgpr), meanpr, grid, cache, gcache, gcachep)
(ftot = 1.1, K = VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}([0.9944957386363635 -0.005326704545454545 … 0.00532670454545453 0.0055042613636363605; -0.005326704545454545 0.9948393969941347 … 0.0051606030058651 0.005326704545454565; … ; 0.00532670454545453 0.0051606030058651 … 0.9948393969941348 -0.005326704545454544; 0.0055042613636363605 0.005326704545454565 … -0.005326704545454544 0.9944957386363638], (32, 32)), meanpr = [-7.55422122563378 -6.823167558636962 … -6.823167558636962 -7.55422122563378; -6.823167558636962 -6.092113891640146 … -6.092113891640146 -6.823167558636962; … ; -6.823167558636962 -6.092113891640146 … -6.092113891640146 -6.823167558636962; -7.55422122563378 -6.823167558636962 … -6.823167558636962 -7.55422122563378], grid = GriddedKeys{(:X, :Y)}
	X: LinRange{Float64}(-3.5224744018114725e-10, 3.5224744018114725e-10, 32)
	Y: LinRange{Float64}(-3.5224744018114725e-10, 3.5224744018114725e-10, 32)
, cache = VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}(VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}(NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}(1, 4, 2.0, :kaiser_bessel, AbstractNFFTs.POLYNOMIAL, true, false, true, 0x00000000), [1.2065151329523807e9 819926.8691406249 … 4.472903950222221e9 -15688.352240668402; -3.8185738483809524e9 -1.7123166484375e6 … -6.2063795222222224e7 119336.55902777778]), NFFTPlan with 284 sampling points for an input array of size(32, 32) and an output array of size(284,) with dims 1:2, ComplexF64[0.9304270810105981 - 0.1755525366839422im, 0.9999999857224657 - 6.37119066797087e-5im, 0.6968566352375541 - 0.5221791050707502im, 0.9304860627318433 + 0.17550202040945986im, 0.8028821861479044 + 0.3950138433443315im, 0.696909133455242 + 0.522149977275377im, 0.9999999859666728 - 6.454986502401626e-5im, 0.690100908057055 - 0.5259755617707915im, 0.9317353097306306 - 0.1646343527093735im, 0.9317937328057632 + 0.1645830409592267im  …  0.8346090489457091 + 0.1795789988845693im, 0.834610800886739 + 0.17958571977595708im, 0.9181402138641586 + 0.09206131734480663im, 0.9120864418553256 + 0.2942607208902332im, 0.9120838428437319 + 0.29426742718619836im, 0.9559925804493471 + 0.09847981391058878im, 0.8881131676530709 + 0.289296040490046im, 0.9924358388500316 + 0.0028007013218525043im, 0.8881114735852544 + 0.28930284794286826im, 0.9999999999233905 + 7.3999334332420275e-6im], BSplinePulse{3}(), [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]), gcache = JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}(sparse([1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [1, 1, 1, 3, 4, 4, 5, 5, 5, 7, 8, 8, 9, 9, 9, 11, 12, 12, 13, 13, 13, 15, 16, 16, 17, 17, 17, 19, 20, 20, 22, 23, 23, 24, 24, 24, 26, 27, 27, 28, 28, 28, 30, 31, 31, 32, 32, 32, 34, 35, 35, 36, 36, 36, 36, 38, 39, 39, 40, 40, 40, 41, 41, 41, 41, 43, 44, 44, 45, 45, 45, 46, 46, 46, 46, 48, 49, 49, 50, 50, 50, 51, 51, 51, 51, 51, 52, 53, 53, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 57, 57, 57, 58, 58, 59, 59, 59, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 64, 64, 64, 64, 64, 64, 65, 65, 66, 66, 66, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 71, 71, 71, 71, 71, 71, 72, 72, 73, 73, 73, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 78, 78, 78, 78, 78, 79, 79, 80, 80, 80, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 86, 86, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 95, 95, 96, 96, 98, 98, 98, 99, 100, 100, 100, 100, 100, 101, 101, 102, 102, 102, 104, 104, 104, 104, 105, 106, 106, 106, 106, 106, 107, 107, 108, 108, 108, 110, 110, 110, 110, 111, 112, 112, 112, 112, 112, 113, 113, 114, 114, 114, 116, 116, 116, 116, 117, 118, 118, 118, 118, 118, 119, 119, 120, 120, 120, 122, 122, 122, 122, 123, 124, 124, 124, 124, 124, 125, 125, 126, 126, 126, 128, 128, 128, 128, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [2, 2, 2, 3, 3, 4, 6, 6, 6, 7, 7, 8, 10, 10, 10, 11, 11, 12, 14, 14, 14, 15, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 22, 25, 25, 25, 26, 26, 27, 29, 29, 29, 30, 30, 31, 33, 33, 33, 34, 34, 35, 37, 37, 37, 37, 38, 38, 38, 39, 39, 40, 42, 42, 42, 42, 43, 43, 43, 44, 44, 45, 47, 47, 47, 47, 48, 48, 48, 49, 49, 50, 52, 52, 52, 52, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 56, 58, 58, 58, 58, 59, 59, 59, 60, 60, 60, 60, 60, 60, 61, 61, 62, 63, 63, 63, 63, 63, 65, 65, 65, 65, 66, 66, 66, 67, 67, 67, 67, 67, 67, 68, 68, 69, 70, 70, 70, 70, 70, 72, 72, 72, 72, 73, 73, 73, 74, 74, 74, 74, 74, 74, 75, 75, 76, 77, 77, 77, 77, 77, 79, 79, 79, 80, 80, 81, 81, 81, 81, 81, 82, 83, 83, 83, 83, 85, 85, 85, 86, 86, 87, 87, 87, 87, 88, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 96, 96, 97, 97, 97, 97, 98, 99, 99, 99, 101, 101, 101, 102, 102, 103, 103, 103, 103, 103, 104, 105, 105, 105, 105, 107, 107, 107, 108, 108, 109, 109, 109, 109, 109, 110, 111, 111, 111, 111, 113, 113, 113, 114, 114, 115, 115, 115, 115, 115, 116, 117, 117, 117, 117, 119, 119, 119, 120, 120, 121, 121, 121, 121, 121, 122, 123, 123, 123, 123, 125, 125, 125, 126, 126, 127, 127, 127, 127, 127, 128, 129, 129, 129, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AA, :AP, :LM, :PV, :AA, :AP, :LM, :PV, :AA, :AP  …  :AZ, :JC, :LM, :SM, :AA, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AA), (0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AA), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AA), (1.516666665673256, :AP)  …  (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AA), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), Fill(NoReference(), 25)), gcachep = JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}(Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([4, 5, 6, 10, 11, 12, 16, 17, 18, 22  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [2, 3, 3, 5, 6, 6, 8, 9, 9, 11, 12, 12, 14, 15, 15, 16, 17, 17, 19, 20, 20, 22, 23, 23, 25, 26, 26, 28, 29, 29, 30, 30, 30, 32, 33, 33, 34, 34, 34, 36, 37, 37, 38, 38, 38, 39, 40, 40, 42, 42, 42, 43, 43, 43, 43, 44, 44, 45, 45, 45, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 50, 50, 51, 51, 51, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 56, 56, 57, 57, 57, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 62, 62, 63, 63, 63, 65, 65, 65, 65, 66, 67, 68, 68, 70, 70, 70, 71, 71, 72, 72, 72, 74, 74, 74, 74, 75, 76, 76, 78, 78, 78, 79, 80, 80, 81, 81, 81, 83, 83, 83, 83, 84, 85, 85, 86, 86, 86, 88, 88, 88, 88, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 96, 96, 96, 98, 98, 98, 98, 99, 100, 100, 101, 101, 101, 103, 103, 103, 103, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [1, 1, 1, 2, 2, 3, 4, 4, 4, 5, 5, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 11, 11, 12, 13, 13, 13, 14, 14, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 21, 22, 22, 23, 24, 24, 24, 25, 25, 26, 27, 27, 27, 27, 28, 28, 28, 29, 29, 30, 31, 31, 31, 31, 32, 32, 32, 33, 33, 34, 35, 35, 35, 35, 36, 36, 36, 37, 37, 38, 39, 39, 39, 39, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 43, 44, 44, 44, 44, 45, 45, 45, 46, 46, 46, 46, 46, 46, 47, 47, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 52, 52, 52, 52, 52, 52, 53, 53, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 58, 58, 58, 58, 58, 58, 59, 59, 60, 61, 61, 61, 61, 61, 62, 62, 62, 63, 63, 64, 64, 64, 64, 64, 65, 66, 66, 66, 66, 67, 67, 67, 68, 68, 69, 69, 69, 69, 70, 71, 71, 71, 72, 72, 73, 73, 73, 73, 73, 74, 75, 75, 75, 75, 76, 76, 77, 77, 77, 77, 78, 79, 79, 79, 80, 80, 80, 81, 81, 82, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 85, 85, 86, 86, 87, 87, 87, 87, 87, 88, 89, 89, 89, 89, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 95, 95, 95, 96, 96, 97, 97, 97, 97, 97, 98, 99, 99, 99, 99, 100, 100, 100, 101, 101, 102, 102, 102, 102, 102, 103, 104, 104, 104, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AP, :LM, :PV, :AP, :LM, :PV, :AP, :LM, :PV, :AP  …  :AP, :AZ, :JC, :LM, :SM, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256, 1.516666665673256, 1.816666603088379  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AP), (1.516666665673256, :LM), (1.516666665673256, :PV), (1.816666603088379, :AP)  …  (7.7166666984558105, :AP), (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), SingleReference{FixedSeg{ComplexF64}}[SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AP, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))  …  SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))]))

We will also fix the total flux to be the observed value 1.1. This is because total flux is degenerate with a global shift in the gain amplitudes making the problem degenerate. To fix this we use the observed total flux as our value.

Moving onto our prior, we first focus on the instrument model priors. Each station requires its own prior on both the amplitudes and phases. For the amplitudes we assume that the gains are apriori well calibrated around unit gains (or 0 log gain amplitudes) which corresponds to no instrument corruption. The gain dispersion is then set to 10% for all stations except LMT, representing that we expect 10% deviations from scan-to-scan. For LMT we let the prior expand to 100% due to the known pointing issues LMT had in 2017.

using Distributions
using DistributionsAD
distamp = station_tuple(dvis, Normal(0.0, 0.1); LM = Normal(1.0))
(AA = Distributions.Normal{Float64}(μ=0.0, σ=0.1), AP = Distributions.Normal{Float64}(μ=0.0, σ=0.1), AZ = Distributions.Normal{Float64}(μ=0.0, σ=0.1), JC = Distributions.Normal{Float64}(μ=0.0, σ=0.1), LM = Distributions.Normal{Float64}(μ=1.0, σ=1.0), PV = Distributions.Normal{Float64}(μ=0.0, σ=0.1), SM = Distributions.Normal{Float64}(μ=0.0, σ=0.1))

For the phases, as mentioned above, we will use a segmented gain prior. This means that rather than the parameters being directly the gains, we fit the first gain for each site, and then the other parameters are the segmented gains compared to the previous time. To model this we break the gain phase prior into two parts. The first is the prior for the first observing timestamp of each site, distphase0, and the second is the prior for segmented gain ϵₜ from time i to i+1, given by distphase. For the EHT, we are dealing with pre-2*rand(rng, ndim) .- 1.5calibrated data, so often, the gain phase jumps from scan to scan are minor. As such, we can put a more informative prior on distphase.

Warning

We use AA (ALMA) as a reference station so we do not have to specify a gain prior for it.

distphase = station_tuple(dvis, DiagonalVonMises(0.0, inv(π^2)))
(AA = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), AP = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), AZ = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), JC = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), LM = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), PV = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688), SM = VLBIImagePriors.DiagonalVonMises{Float64, Float64, Float64}(μ=0.0, κ=0.10132118364233778, lnorm=-1.739120733481688))

In addition we want a reasonable guess for what the resolution of our image should be. For radio astronomy this is given by roughly the longest baseline in the image. To put this into pixel space we then divide by the pixel size.

beam = beamsize(dvis)
rat = (beam/(step(grid.X)))
5.279832635689346

To make the Gaussian Markov random field efficient we first precompute a bunch of quantities that allow us to scale things linearly with the number of image pixels. This drastically improves the usual N^3 scaling you get from usual Gaussian Processes.

crcache = MarkovRandomFieldCache(meanpr)
VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}(sparse([1, 2, 32, 33, 993, 1, 2, 3, 34, 994  …  31, 991, 1022, 1023, 1024, 32, 992, 993, 1023, 1024], [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 19, 19, 19, 19, 19, 20, 20, 20, 20, 20, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 24, 24, 24, 24, 24, 25, 25, 25, 25, 25, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, 30, 30, 30, 30, 30, 31, 31, 31, 31, 31, 32, 32, 32, 32, 32, 33, 33, 33, 33, 33, 34, 34, 34, 34, 34, 35, 35, 35, 35, 35, 36, 36, 36, 36, 36, 37, 37, 37, 37, 37, 38, 38, 38, 38, 38, 39, 39, 39, 39, 39, 40, 40, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 42, 42, 42, 43, 43, 43, 43, 43, 44, 44, 44, 44, 44, 45, 45, 45, 45, 45, 46, 46, 46, 46, 46, 47, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 50, 51, 51, 51, 51, 51, 52, 52, 52, 52, 52, 53, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 56, 57, 57, 57, 57, 57, 58, 58, 58, 58, 58, 59, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 63, 63, 63, 63, 64, 64, 64, 64, 64, 65, 65, 65, 65, 65, 66, 66, 66, 66, 66, 67, 67, 67, 67, 67, 68, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 70, 70, 70, 70, 71, 71, 71, 71, 71, 72, 72, 72, 72, 72, 73, 73, 73, 73, 73, 74, 74, 74, 74, 74, 75, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 77, 77, 77, 77, 78, 78, 78, 78, 78, 79, 79, 79, 79, 79, 80, 80, 80, 80, 80, 81, 81, 81, 81, 81, 82, 82, 82, 82, 82, 83, 83, 83, 83, 83, 84, 84, 84, 84, 84, 85, 85, 85, 85, 85, 86, 86, 86, 86, 86, 87, 87, 87, 87, 87, 88, 88, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 90, 90, 90, 91, 91, 91, 91, 91, 92, 92, 92, 92, 92, 93, 93, 93, 93, 93, 94, 94, 94, 94, 94, 95, 95, 95, 95, 95, 96, 96, 96, 96, 96, 97, 97, 97, 97, 97, 98, 98, 98, 98, 98, 99, 99, 99, 99, 99, 100, 100, 100, 100, 100, 101, 101, 101, 101, 101, 102, 102, 102, 102, 102, 103, 103, 103, 103, 103, 104, 104, 104, 104, 104, 105, 105, 105, 105, 105, 106, 106, 106, 106, 106, 107, 107, 107, 107, 107, 108, 108, 108, 108, 108, 109, 109, 109, 109, 109, 110, 110, 110, 110, 110, 111, 111, 111, 111, 111, 112, 112, 112, 112, 112, 113, 113, 113, 113, 113, 114, 114, 114, 114, 114, 115, 115, 115, 115, 115, 116, 116, 116, 116, 116, 117, 117, 117, 117, 117, 118, 118, 118, 118, 118, 119, 119, 119, 119, 119, 120, 120, 120, 120, 120, 121, 121, 121, 121, 121, 122, 122, 122, 122, 122, 123, 123, 123, 123, 123, 124, 124, 124, 124, 124, 125, 125, 125, 125, 125, 126, 126, 126, 126, 126, 127, 127, 127, 127, 127, 128, 128, 128, 128, 128, 129, 129, 129, 129, 129, 130, 130, 130, 130, 130, 131, 131, 131, 131, 131, 132, 132, 132, 132, 132, 133, 133, 133, 133, 133, 134, 134, 134, 134, 134, 135, 135, 135, 135, 135, 136, 136, 136, 136, 136, 137, 137, 137, 137, 137, 138, 138, 138, 138, 138, 139, 139, 139, 139, 139, 140, 140, 140, 140, 140, 141, 141, 141, 141, 141, 142, 142, 142, 142, 142, 143, 143, 143, 143, 143, 144, 144, 144, 144, 144, 145, 145, 145, 145, 145, 146, 146, 146, 146, 146, 147, 147, 147, 147, 147, 148, 148, 148, 148, 148, 149, 149, 149, 149, 149, 150, 150, 150, 150, 150, 151, 151, 151, 151, 151, 152, 152, 152, 152, 152, 153, 153, 153, 153, 153, 154, 154, 154, 154, 154, 155, 155, 155, 155, 155, 156, 156, 156, 156, 156, 157, 157, 157, 157, 157, 158, 158, 158, 158, 158, 159, 159, 159, 159, 159, 160, 160, 160, 160, 160, 161, 161, 161, 161, 161, 162, 162, 162, 162, 162, 163, 163, 163, 163, 163, 164, 164, 164, 164, 164, 165, 165, 165, 165, 165, 166, 166, 166, 166, 166, 167, 167, 167, 167, 167, 168, 168, 168, 168, 168, 169, 169, 169, 169, 169, 170, 170, 170, 170, 170, 171, 171, 171, 171, 171, 172, 172, 172, 172, 172, 173, 173, 173, 173, 173, 174, 174, 174, 174, 174, 175, 175, 175, 175, 175, 176, 176, 176, 176, 176, 177, 177, 177, 177, 177, 178, 178, 178, 178, 178, 179, 179, 179, 179, 179, 180, 180, 180, 180, 180, 181, 181, 181, 181, 181, 182, 182, 182, 182, 182, 183, 183, 183, 183, 183, 184, 184, 184, 184, 184, 185, 185, 185, 185, 185, 186, 186, 186, 186, 186, 187, 187, 187, 187, 187, 188, 188, 188, 188, 188, 189, 189, 189, 189, 189, 190, 190, 190, 190, 190, 191, 191, 191, 191, 191, 192, 192, 192, 192, 192, 193, 193, 193, 193, 193, 194, 194, 194, 194, 194, 195, 195, 195, 195, 195, 196, 196, 196, 196, 196, 197, 197, 197, 197, 197, 198, 198, 198, 198, 198, 199, 199, 199, 199, 199, 200, 200, 200, 200, 200, 201, 201, 201, 201, 201, 202, 202, 202, 202, 202, 203, 203, 203, 203, 203, 204, 204, 204, 204, 204, 205, 205, 205, 205, 205, 206, 206, 206, 206, 206, 207, 207, 207, 207, 207, 208, 208, 208, 208, 208, 209, 209, 209, 209, 209, 210, 210, 210, 210, 210, 211, 211, 211, 211, 211, 212, 212, 212, 212, 212, 213, 213, 213, 213, 213, 214, 214, 214, 214, 214, 215, 215, 215, 215, 215, 216, 216, 216, 216, 216, 217, 217, 217, 217, 217, 218, 218, 218, 218, 218, 219, 219, 219, 219, 219, 220, 220, 220, 220, 220, 221, 221, 221, 221, 221, 222, 222, 222, 222, 222, 223, 223, 223, 223, 223, 224, 224, 224, 224, 224, 225, 225, 225, 225, 225, 226, 226, 226, 226, 226, 227, 227, 227, 227, 227, 228, 228, 228, 228, 228, 229, 229, 229, 229, 229, 230, 230, 230, 230, 230, 231, 231, 231, 231, 231, 232, 232, 232, 232, 232, 233, 233, 233, 233, 233, 234, 234, 234, 234, 234, 235, 235, 235, 235, 235, 236, 236, 236, 236, 236, 237, 237, 237, 237, 237, 238, 238, 238, 238, 238, 239, 239, 239, 239, 239, 240, 240, 240, 240, 240, 241, 241, 241, 241, 241, 242, 242, 242, 242, 242, 243, 243, 243, 243, 243, 244, 244, 244, 244, 244, 245, 245, 245, 245, 245, 246, 246, 246, 246, 246, 247, 247, 247, 247, 247, 248, 248, 248, 248, 248, 249, 249, 249, 249, 249, 250, 250, 250, 250, 250, 251, 251, 251, 251, 251, 252, 252, 252, 252, 252, 253, 253, 253, 253, 253, 254, 254, 254, 254, 254, 255, 255, 255, 255, 255, 256, 256, 256, 256, 256, 257, 257, 257, 257, 257, 258, 258, 258, 258, 258, 259, 259, 259, 259, 259, 260, 260, 260, 260, 260, 261, 261, 261, 261, 261, 262, 262, 262, 262, 262, 263, 263, 263, 263, 263, 264, 264, 264, 264, 264, 265, 265, 265, 265, 265, 266, 266, 266, 266, 266, 267, 267, 267, 267, 267, 268, 268, 268, 268, 268, 269, 269, 269, 269, 269, 270, 270, 270, 270, 270, 271, 271, 271, 271, 271, 272, 272, 272, 272, 272, 273, 273, 273, 273, 273, 274, 274, 274, 274, 274, 275, 275, 275, 275, 275, 276, 276, 276, 276, 276, 277, 277, 277, 277, 277, 278, 278, 278, 278, 278, 279, 279, 279, 279, 279, 280, 280, 280, 280, 280, 281, 281, 281, 281, 281, 282, 282, 282, 282, 282, 283, 283, 283, 283, 283, 284, 284, 284, 284, 284, 285, 285, 285, 285, 285, 286, 286, 286, 286, 286, 287, 287, 287, 287, 287, 288, 288, 288, 288, 288, 289, 289, 289, 289, 289, 290, 290, 290, 290, 290, 291, 291, 291, 291, 291, 292, 292, 292, 292, 292, 293, 293, 293, 293, 293, 294, 294, 294, 294, 294, 295, 295, 295, 295, 295, 296, 296, 296, 296, 296, 297, 297, 297, 297, 297, 298, 298, 298, 298, 298, 299, 299, 299, 299, 299, 300, 300, 300, 300, 300, 301, 301, 301, 301, 301, 302, 302, 302, 302, 302, 303, 303, 303, 303, 303, 304, 304, 304, 304, 304, 305, 305, 305, 305, 305, 306, 306, 306, 306, 306, 307, 307, 307, 307, 307, 308, 308, 308, 308, 308, 309, 309, 309, 309, 309, 310, 310, 310, 310, 310, 311, 311, 311, 311, 311, 312, 312, 312, 312, 312, 313, 313, 313, 313, 313, 314, 314, 314, 314, 314, 315, 315, 315, 315, 315, 316, 316, 316, 316, 316, 317, 317, 317, 317, 317, 318, 318, 318, 318, 318, 319, 319, 319, 319, 319, 320, 320, 320, 320, 320, 321, 321, 321, 321, 321, 322, 322, 322, 322, 322, 323, 323, 323, 323, 323, 324, 324, 324, 324, 324, 325, 325, 325, 325, 325, 326, 326, 326, 326, 326, 327, 327, 327, 327, 327, 328, 328, 328, 328, 328, 329, 329, 329, 329, 329, 330, 330, 330, 330, 330, 331, 331, 331, 331, 331, 332, 332, 332, 332, 332, 333, 333, 333, 333, 333, 334, 334, 334, 334, 334, 335, 335, 335, 335, 335, 336, 336, 336, 336, 336, 337, 337, 337, 337, 337, 338, 338, 338, 338, 338, 339, 339, 339, 339, 339, 340, 340, 340, 340, 340, 341, 341, 341, 341, 341, 342, 342, 342, 342, 342, 343, 343, 343, 343, 343, 344, 344, 344, 344, 344, 345, 345, 345, 345, 345, 346, 346, 346, 346, 346, 347, 347, 347, 347, 347, 348, 348, 348, 348, 348, 349, 349, 349, 349, 349, 350, 350, 350, 350, 350, 351, 351, 351, 351, 351, 352, 352, 352, 352, 352, 353, 353, 353, 353, 353, 354, 354, 354, 354, 354, 355, 355, 355, 355, 355, 356, 356, 356, 356, 356, 357, 357, 357, 357, 357, 358, 358, 358, 358, 358, 359, 359, 359, 359, 359, 360, 360, 360, 360, 360, 361, 361, 361, 361, 361, 362, 362, 362, 362, 362, 363, 363, 363, 363, 363, 364, 364, 364, 364, 364, 365, 365, 365, 365, 365, 366, 366, 366, 366, 366, 367, 367, 367, 367, 367, 368, 368, 368, 368, 368, 369, 369, 369, 369, 369, 370, 370, 370, 370, 370, 371, 371, 371, 371, 371, 372, 372, 372, 372, 372, 373, 373, 373, 373, 373, 374, 374, 374, 374, 374, 375, 375, 375, 375, 375, 376, 376, 376, 376, 376, 377, 377, 377, 377, 377, 378, 378, 378, 378, 378, 379, 379, 379, 379, 379, 380, 380, 380, 380, 380, 381, 381, 381, 381, 381, 382, 382, 382, 382, 382, 383, 383, 383, 383, 383, 384, 384, 384, 384, 384, 385, 385, 385, 385, 385, 386, 386, 386, 386, 386, 387, 387, 387, 387, 387, 388, 388, 388, 388, 388, 389, 389, 389, 389, 389, 390, 390, 390, 390, 390, 391, 391, 391, 391, 391, 392, 392, 392, 392, 392, 393, 393, 393, 393, 393, 394, 394, 394, 394, 394, 395, 395, 395, 395, 395, 396, 396, 396, 396, 396, 397, 397, 397, 397, 397, 398, 398, 398, 398, 398, 399, 399, 399, 399, 399, 400, 400, 400, 400, 400, 401, 401, 401, 401, 401, 402, 402, 402, 402, 402, 403, 403, 403, 403, 403, 404, 404, 404, 404, 404, 405, 405, 405, 405, 405, 406, 406, 406, 406, 406, 407, 407, 407, 407, 407, 408, 408, 408, 408, 408, 409, 409, 409, 409, 409, 410, 410, 410, 410, 410, 411, 411, 411, 411, 411, 412, 412, 412, 412, 412, 413, 413, 413, 413, 413, 414, 414, 414, 414, 414, 415, 415, 415, 415, 415, 416, 416, 416, 416, 416, 417, 417, 417, 417, 417, 418, 418, 418, 418, 418, 419, 419, 419, 419, 419, 420, 420, 420, 420, 420, 421, 421, 421, 421, 421, 422, 422, 422, 422, 422, 423, 423, 423, 423, 423, 424, 424, 424, 424, 424, 425, 425, 425, 425, 425, 426, 426, 426, 426, 426, 427, 427, 427, 427, 427, 428, 428, 428, 428, 428, 429, 429, 429, 429, 429, 430, 430, 430, 430, 430, 431, 431, 431, 431, 431, 432, 432, 432, 432, 432, 433, 433, 433, 433, 433, 434, 434, 434, 434, 434, 435, 435, 435, 435, 435, 436, 436, 436, 436, 436, 437, 437, 437, 437, 437, 438, 438, 438, 438, 438, 439, 439, 439, 439, 439, 440, 440, 440, 440, 440, 441, 441, 441, 441, 441, 442, 442, 442, 442, 442, 443, 443, 443, 443, 443, 444, 444, 444, 444, 444, 445, 445, 445, 445, 445, 446, 446, 446, 446, 446, 447, 447, 447, 447, 447, 448, 448, 448, 448, 448, 449, 449, 449, 449, 449, 450, 450, 450, 450, 450, 451, 451, 451, 451, 451, 452, 452, 452, 452, 452, 453, 453, 453, 453, 453, 454, 454, 454, 454, 454, 455, 455, 455, 455, 455, 456, 456, 456, 456, 456, 457, 457, 457, 457, 457, 458, 458, 458, 458, 458, 459, 459, 459, 459, 459, 460, 460, 460, 460, 460, 461, 461, 461, 461, 461, 462, 462, 462, 462, 462, 463, 463, 463, 463, 463, 464, 464, 464, 464, 464, 465, 465, 465, 465, 465, 466, 466, 466, 466, 466, 467, 467, 467, 467, 467, 468, 468, 468, 468, 468, 469, 469, 469, 469, 469, 470, 470, 470, 470, 470, 471, 471, 471, 471, 471, 472, 472, 472, 472, 472, 473, 473, 473, 473, 473, 474, 474, 474, 474, 474, 475, 475, 475, 475, 475, 476, 476, 476, 476, 476, 477, 477, 477, 477, 477, 478, 478, 478, 478, 478, 479, 479, 479, 479, 479, 480, 480, 480, 480, 480, 481, 481, 481, 481, 481, 482, 482, 482, 482, 482, 483, 483, 483, 483, 483, 484, 484, 484, 484, 484, 485, 485, 485, 485, 485, 486, 486, 486, 486, 486, 487, 487, 487, 487, 487, 488, 488, 488, 488, 488, 489, 489, 489, 489, 489, 490, 490, 490, 490, 490, 491, 491, 491, 491, 491, 492, 492, 492, 492, 492, 493, 493, 493, 493, 493, 494, 494, 494, 494, 494, 495, 495, 495, 495, 495, 496, 496, 496, 496, 496, 497, 497, 497, 497, 497, 498, 498, 498, 498, 498, 499, 499, 499, 499, 499, 500, 500, 500, 500, 500, 501, 501, 501, 501, 501, 502, 502, 502, 502, 502, 503, 503, 503, 503, 503, 504, 504, 504, 504, 504, 505, 505, 505, 505, 505, 506, 506, 506, 506, 506, 507, 507, 507, 507, 507, 508, 508, 508, 508, 508, 509, 509, 509, 509, 509, 510, 510, 510, 510, 510, 511, 511, 511, 511, 511, 512, 512, 512, 512, 512, 513, 513, 513, 513, 513, 514, 514, 514, 514, 514, 515, 515, 515, 515, 515, 516, 516, 516, 516, 516, 517, 517, 517, 517, 517, 518, 518, 518, 518, 518, 519, 519, 519, 519, 519, 520, 520, 520, 520, 520, 521, 521, 521, 521, 521, 522, 522, 522, 522, 522, 523, 523, 523, 523, 523, 524, 524, 524, 524, 524, 525, 525, 525, 525, 525, 526, 526, 526, 526, 526, 527, 527, 527, 527, 527, 528, 528, 528, 528, 528, 529, 529, 529, 529, 529, 530, 530, 530, 530, 530, 531, 531, 531, 531, 531, 532, 532, 532, 532, 532, 533, 533, 533, 533, 533, 534, 534, 534, 534, 534, 535, 535, 535, 535, 535, 536, 536, 536, 536, 536, 537, 537, 537, 537, 537, 538, 538, 538, 538, 538, 539, 539, 539, 539, 539, 540, 540, 540, 540, 540, 541, 541, 541, 541, 541, 542, 542, 542, 542, 542, 543, 543, 543, 543, 543, 544, 544, 544, 544, 544, 545, 545, 545, 545, 545, 546, 546, 546, 546, 546, 547, 547, 547, 547, 547, 548, 548, 548, 548, 548, 549, 549, 549, 549, 549, 550, 550, 550, 550, 550, 551, 551, 551, 551, 551, 552, 552, 552, 552, 552, 553, 553, 553, 553, 553, 554, 554, 554, 554, 554, 555, 555, 555, 555, 555, 556, 556, 556, 556, 556, 557, 557, 557, 557, 557, 558, 558, 558, 558, 558, 559, 559, 559, 559, 559, 560, 560, 560, 560, 560, 561, 561, 561, 561, 561, 562, 562, 562, 562, 562, 563, 563, 563, 563, 563, 564, 564, 564, 564, 564, 565, 565, 565, 565, 565, 566, 566, 566, 566, 566, 567, 567, 567, 567, 567, 568, 568, 568, 568, 568, 569, 569, 569, 569, 569, 570, 570, 570, 570, 570, 571, 571, 571, 571, 571, 572, 572, 572, 572, 572, 573, 573, 573, 573, 573, 574, 574, 574, 574, 574, 575, 575, 575, 575, 575, 576, 576, 576, 576, 576, 577, 577, 577, 577, 577, 578, 578, 578, 578, 578, 579, 579, 579, 579, 579, 580, 580, 580, 580, 580, 581, 581, 581, 581, 581, 582, 582, 582, 582, 582, 583, 583, 583, 583, 583, 584, 584, 584, 584, 584, 585, 585, 585, 585, 585, 586, 586, 586, 586, 586, 587, 587, 587, 587, 587, 588, 588, 588, 588, 588, 589, 589, 589, 589, 589, 590, 590, 590, 590, 590, 591, 591, 591, 591, 591, 592, 592, 592, 592, 592, 593, 593, 593, 593, 593, 594, 594, 594, 594, 594, 595, 595, 595, 595, 595, 596, 596, 596, 596, 596, 597, 597, 597, 597, 597, 598, 598, 598, 598, 598, 599, 599, 599, 599, 599, 600, 600, 600, 600, 600, 601, 601, 601, 601, 601, 602, 602, 602, 602, 602, 603, 603, 603, 603, 603, 604, 604, 604, 604, 604, 605, 605, 605, 605, 605, 606, 606, 606, 606, 606, 607, 607, 607, 607, 607, 608, 608, 608, 608, 608, 609, 609, 609, 609, 609, 610, 610, 610, 610, 610, 611, 611, 611, 611, 611, 612, 612, 612, 612, 612, 613, 613, 613, 613, 613, 614, 614, 614, 614, 614, 615, 615, 615, 615, 615, 616, 616, 616, 616, 616, 617, 617, 617, 617, 617, 618, 618, 618, 618, 618, 619, 619, 619, 619, 619, 620, 620, 620, 620, 620, 621, 621, 621, 621, 621, 622, 622, 622, 622, 622, 623, 623, 623, 623, 623, 624, 624, 624, 624, 624, 625, 625, 625, 625, 625, 626, 626, 626, 626, 626, 627, 627, 627, 627, 627, 628, 628, 628, 628, 628, 629, 629, 629, 629, 629, 630, 630, 630, 630, 630, 631, 631, 631, 631, 631, 632, 632, 632, 632, 632, 633, 633, 633, 633, 633, 634, 634, 634, 634, 634, 635, 635, 635, 635, 635, 636, 636, 636, 636, 636, 637, 637, 637, 637, 637, 638, 638, 638, 638, 638, 639, 639, 639, 639, 639, 640, 640, 640, 640, 640, 641, 641, 641, 641, 641, 642, 642, 642, 642, 642, 643, 643, 643, 643, 643, 644, 644, 644, 644, 644, 645, 645, 645, 645, 645, 646, 646, 646, 646, 646, 647, 647, 647, 647, 647, 648, 648, 648, 648, 648, 649, 649, 649, 649, 649, 650, 650, 650, 650, 650, 651, 651, 651, 651, 651, 652, 652, 652, 652, 652, 653, 653, 653, 653, 653, 654, 654, 654, 654, 654, 655, 655, 655, 655, 655, 656, 656, 656, 656, 656, 657, 657, 657, 657, 657, 658, 658, 658, 658, 658, 659, 659, 659, 659, 659, 660, 660, 660, 660, 660, 661, 661, 661, 661, 661, 662, 662, 662, 662, 662, 663, 663, 663, 663, 663, 664, 664, 664, 664, 664, 665, 665, 665, 665, 665, 666, 666, 666, 666, 666, 667, 667, 667, 667, 667, 668, 668, 668, 668, 668, 669, 669, 669, 669, 669, 670, 670, 670, 670, 670, 671, 671, 671, 671, 671, 672, 672, 672, 672, 672, 673, 673, 673, 673, 673, 674, 674, 674, 674, 674, 675, 675, 675, 675, 675, 676, 676, 676, 676, 676, 677, 677, 677, 677, 677, 678, 678, 678, 678, 678, 679, 679, 679, 679, 679, 680, 680, 680, 680, 680, 681, 681, 681, 681, 681, 682, 682, 682, 682, 682, 683, 683, 683, 683, 683, 684, 684, 684, 684, 684, 685, 685, 685, 685, 685, 686, 686, 686, 686, 686, 687, 687, 687, 687, 687, 688, 688, 688, 688, 688, 689, 689, 689, 689, 689, 690, 690, 690, 690, 690, 691, 691, 691, 691, 691, 692, 692, 692, 692, 692, 693, 693, 693, 693, 693, 694, 694, 694, 694, 694, 695, 695, 695, 695, 695, 696, 696, 696, 696, 696, 697, 697, 697, 697, 697, 698, 698, 698, 698, 698, 699, 699, 699, 699, 699, 700, 700, 700, 700, 700, 701, 701, 701, 701, 701, 702, 702, 702, 702, 702, 703, 703, 703, 703, 703, 704, 704, 704, 704, 704, 705, 705, 705, 705, 705, 706, 706, 706, 706, 706, 707, 707, 707, 707, 707, 708, 708, 708, 708, 708, 709, 709, 709, 709, 709, 710, 710, 710, 710, 710, 711, 711, 711, 711, 711, 712, 712, 712, 712, 712, 713, 713, 713, 713, 713, 714, 714, 714, 714, 714, 715, 715, 715, 715, 715, 716, 716, 716, 716, 716, 717, 717, 717, 717, 717, 718, 718, 718, 718, 718, 719, 719, 719, 719, 719, 720, 720, 720, 720, 720, 721, 721, 721, 721, 721, 722, 722, 722, 722, 722, 723, 723, 723, 723, 723, 724, 724, 724, 724, 724, 725, 725, 725, 725, 725, 726, 726, 726, 726, 726, 727, 727, 727, 727, 727, 728, 728, 728, 728, 728, 729, 729, 729, 729, 729, 730, 730, 730, 730, 730, 731, 731, 731, 731, 731, 732, 732, 732, 732, 732, 733, 733, 733, 733, 733, 734, 734, 734, 734, 734, 735, 735, 735, 735, 735, 736, 736, 736, 736, 736, 737, 737, 737, 737, 737, 738, 738, 738, 738, 738, 739, 739, 739, 739, 739, 740, 740, 740, 740, 740, 741, 741, 741, 741, 741, 742, 742, 742, 742, 742, 743, 743, 743, 743, 743, 744, 744, 744, 744, 744, 745, 745, 745, 745, 745, 746, 746, 746, 746, 746, 747, 747, 747, 747, 747, 748, 748, 748, 748, 748, 749, 749, 749, 749, 749, 750, 750, 750, 750, 750, 751, 751, 751, 751, 751, 752, 752, 752, 752, 752, 753, 753, 753, 753, 753, 754, 754, 754, 754, 754, 755, 755, 755, 755, 755, 756, 756, 756, 756, 756, 757, 757, 757, 757, 757, 758, 758, 758, 758, 758, 759, 759, 759, 759, 759, 760, 760, 760, 760, 760, 761, 761, 761, 761, 761, 762, 762, 762, 762, 762, 763, 763, 763, 763, 763, 764, 764, 764, 764, 764, 765, 765, 765, 765, 765, 766, 766, 766, 766, 766, 767, 767, 767, 767, 767, 768, 768, 768, 768, 768, 769, 769, 769, 769, 769, 770, 770, 770, 770, 770, 771, 771, 771, 771, 771, 772, 772, 772, 772, 772, 773, 773, 773, 773, 773, 774, 774, 774, 774, 774, 775, 775, 775, 775, 775, 776, 776, 776, 776, 776, 777, 777, 777, 777, 777, 778, 778, 778, 778, 778, 779, 779, 779, 779, 779, 780, 780, 780, 780, 780, 781, 781, 781, 781, 781, 782, 782, 782, 782, 782, 783, 783, 783, 783, 783, 784, 784, 784, 784, 784, 785, 785, 785, 785, 785, 786, 786, 786, 786, 786, 787, 787, 787, 787, 787, 788, 788, 788, 788, 788, 789, 789, 789, 789, 789, 790, 790, 790, 790, 790, 791, 791, 791, 791, 791, 792, 792, 792, 792, 792, 793, 793, 793, 793, 793, 794, 794, 794, 794, 794, 795, 795, 795, 795, 795, 796, 796, 796, 796, 796, 797, 797, 797, 797, 797, 798, 798, 798, 798, 798, 799, 799, 799, 799, 799, 800, 800, 800, 800, 800, 801, 801, 801, 801, 801, 802, 802, 802, 802, 802, 803, 803, 803, 803, 803, 804, 804, 804, 804, 804, 805, 805, 805, 805, 805, 806, 806, 806, 806, 806, 807, 807, 807, 807, 807, 808, 808, 808, 808, 808, 809, 809, 809, 809, 809, 810, 810, 810, 810, 810, 811, 811, 811, 811, 811, 812, 812, 812, 812, 812, 813, 813, 813, 813, 813, 814, 814, 814, 814, 814, 815, 815, 815, 815, 815, 816, 816, 816, 816, 816, 817, 817, 817, 817, 817, 818, 818, 818, 818, 818, 819, 819, 819, 819, 819, 820, 820, 820, 820, 820, 821, 821, 821, 821, 821, 822, 822, 822, 822, 822, 823, 823, 823, 823, 823, 824, 824, 824, 824, 824, 825, 825, 825, 825, 825, 826, 826, 826, 826, 826, 827, 827, 827, 827, 827, 828, 828, 828, 828, 828, 829, 829, 829, 829, 829, 830, 830, 830, 830, 830, 831, 831, 831, 831, 831, 832, 832, 832, 832, 832, 833, 833, 833, 833, 833, 834, 834, 834, 834, 834, 835, 835, 835, 835, 835, 836, 836, 836, 836, 836, 837, 837, 837, 837, 837, 838, 838, 838, 838, 838, 839, 839, 839, 839, 839, 840, 840, 840, 840, 840, 841, 841, 841, 841, 841, 842, 842, 842, 842, 842, 843, 843, 843, 843, 843, 844, 844, 844, 844, 844, 845, 845, 845, 845, 845, 846, 846, 846, 846, 846, 847, 847, 847, 847, 847, 848, 848, 848, 848, 848, 849, 849, 849, 849, 849, 850, 850, 850, 850, 850, 851, 851, 851, 851, 851, 852, 852, 852, 852, 852, 853, 853, 853, 853, 853, 854, 854, 854, 854, 854, 855, 855, 855, 855, 855, 856, 856, 856, 856, 856, 857, 857, 857, 857, 857, 858, 858, 858, 858, 858, 859, 859, 859, 859, 859, 860, 860, 860, 860, 860, 861, 861, 861, 861, 861, 862, 862, 862, 862, 862, 863, 863, 863, 863, 863, 864, 864, 864, 864, 864, 865, 865, 865, 865, 865, 866, 866, 866, 866, 866, 867, 867, 867, 867, 867, 868, 868, 868, 868, 868, 869, 869, 869, 869, 869, 870, 870, 870, 870, 870, 871, 871, 871, 871, 871, 872, 872, 872, 872, 872, 873, 873, 873, 873, 873, 874, 874, 874, 874, 874, 875, 875, 875, 875, 875, 876, 876, 876, 876, 876, 877, 877, 877, 877, 877, 878, 878, 878, 878, 878, 879, 879, 879, 879, 879, 880, 880, 880, 880, 880, 881, 881, 881, 881, 881, 882, 882, 882, 882, 882, 883, 883, 883, 883, 883, 884, 884, 884, 884, 884, 885, 885, 885, 885, 885, 886, 886, 886, 886, 886, 887, 887, 887, 887, 887, 888, 888, 888, 888, 888, 889, 889, 889, 889, 889, 890, 890, 890, 890, 890, 891, 891, 891, 891, 891, 892, 892, 892, 892, 892, 893, 893, 893, 893, 893, 894, 894, 894, 894, 894, 895, 895, 895, 895, 895, 896, 896, 896, 896, 896, 897, 897, 897, 897, 897, 898, 898, 898, 898, 898, 899, 899, 899, 899, 899, 900, 900, 900, 900, 900, 901, 901, 901, 901, 901, 902, 902, 902, 902, 902, 903, 903, 903, 903, 903, 904, 904, 904, 904, 904, 905, 905, 905, 905, 905, 906, 906, 906, 906, 906, 907, 907, 907, 907, 907, 908, 908, 908, 908, 908, 909, 909, 909, 909, 909, 910, 910, 910, 910, 910, 911, 911, 911, 911, 911, 912, 912, 912, 912, 912, 913, 913, 913, 913, 913, 914, 914, 914, 914, 914, 915, 915, 915, 915, 915, 916, 916, 916, 916, 916, 917, 917, 917, 917, 917, 918, 918, 918, 918, 918, 919, 919, 919, 919, 919, 920, 920, 920, 920, 920, 921, 921, 921, 921, 921, 922, 922, 922, 922, 922, 923, 923, 923, 923, 923, 924, 924, 924, 924, 924, 925, 925, 925, 925, 925, 926, 926, 926, 926, 926, 927, 927, 927, 927, 927, 928, 928, 928, 928, 928, 929, 929, 929, 929, 929, 930, 930, 930, 930, 930, 931, 931, 931, 931, 931, 932, 932, 932, 932, 932, 933, 933, 933, 933, 933, 934, 934, 934, 934, 934, 935, 935, 935, 935, 935, 936, 936, 936, 936, 936, 937, 937, 937, 937, 937, 938, 938, 938, 938, 938, 939, 939, 939, 939, 939, 940, 940, 940, 940, 940, 941, 941, 941, 941, 941, 942, 942, 942, 942, 942, 943, 943, 943, 943, 943, 944, 944, 944, 944, 944, 945, 945, 945, 945, 945, 946, 946, 946, 946, 946, 947, 947, 947, 947, 947, 948, 948, 948, 948, 948, 949, 949, 949, 949, 949, 950, 950, 950, 950, 950, 951, 951, 951, 951, 951, 952, 952, 952, 952, 952, 953, 953, 953, 953, 953, 954, 954, 954, 954, 954, 955, 955, 955, 955, 955, 956, 956, 956, 956, 956, 957, 957, 957, 957, 957, 958, 958, 958, 958, 958, 959, 959, 959, 959, 959, 960, 960, 960, 960, 960, 961, 961, 961, 961, 961, 962, 962, 962, 962, 962, 963, 963, 963, 963, 963, 964, 964, 964, 964, 964, 965, 965, 965, 965, 965, 966, 966, 966, 966, 966, 967, 967, 967, 967, 967, 968, 968, 968, 968, 968, 969, 969, 969, 969, 969, 970, 970, 970, 970, 970, 971, 971, 971, 971, 971, 972, 972, 972, 972, 972, 973, 973, 973, 973, 973, 974, 974, 974, 974, 974, 975, 975, 975, 975, 975, 976, 976, 976, 976, 976, 977, 977, 977, 977, 977, 978, 978, 978, 978, 978, 979, 979, 979, 979, 979, 980, 980, 980, 980, 980, 981, 981, 981, 981, 981, 982, 982, 982, 982, 982, 983, 983, 983, 983, 983, 984, 984, 984, 984, 984, 985, 985, 985, 985, 985, 986, 986, 986, 986, 986, 987, 987, 987, 987, 987, 988, 988, 988, 988, 988, 989, 989, 989, 989, 989, 990, 990, 990, 990, 990, 991, 991, 991, 991, 991, 992, 992, 992, 992, 992, 993, 993, 993, 993, 993, 994, 994, 994, 994, 994, 995, 995, 995, 995, 995, 996, 996, 996, 996, 996, 997, 997, 997, 997, 997, 998, 998, 998, 998, 998, 999, 999, 999, 999, 999, 1000, 1000, 1000, 1000, 1000, 1001, 1001, 1001, 1001, 1001, 1002, 1002, 1002, 1002, 1002, 1003, 1003, 1003, 1003, 1003, 1004, 1004, 1004, 1004, 1004, 1005, 1005, 1005, 1005, 1005, 1006, 1006, 1006, 1006, 1006, 1007, 1007, 1007, 1007, 1007, 1008, 1008, 1008, 1008, 1008, 1009, 1009, 1009, 1009, 1009, 1010, 1010, 1010, 1010, 1010, 1011, 1011, 1011, 1011, 1011, 1012, 1012, 1012, 1012, 1012, 1013, 1013, 1013, 1013, 1013, 1014, 1014, 1014, 1014, 1014, 1015, 1015, 1015, 1015, 1015, 1016, 1016, 1016, 1016, 1016, 1017, 1017, 1017, 1017, 1017, 1018, 1018, 1018, 1018, 1018, 1019, 1019, 1019, 1019, 1019, 1020, 1020, 1020, 1020, 1020, 1021, 1021, 1021, 1021, 1021, 1022, 1022, 1022, 1022, 1022, 1023, 1023, 1023, 1023, 1023, 1024, 1024, 1024, 1024, 1024], [4.0, -1.0, -1.0, -1.0, -1.0, -1.0, 4.0, -1.0, -1.0, -1.0  …  -1.0, -1.0, -1.0, 4.0, -1.0, -1.0, -1.0, -1.0, -1.0, 4.0], 1024, 1024), [1.0 0.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 1.0], [0.0 0.03842943919353914 … 0.15224093497742697 0.03842943919353936; 0.03842943919353914 0.07685887838707828 … 0.1906703741709661 0.0768588783870785; … ; 0.15224093497742697 0.1906703741709661 … 0.30448186995485393 0.19067037417096633; 0.03842943919353914 0.07685887838707828 … 0.1906703741709661 0.0768588783870785])

One of the benefits of the Bayesian approach is that we can fit for the hyperparameters of our prior/regularizers unlike traditional RML appraoches. To construct this heirarchical prior we will first make a map that takes in our regularizer hyperparameters and returns the image prior given those hyperparameters.

fmap = let crcache=crcache
    x->GaussMarkovRandomField(x+1, 1.0, crcache)
end
#1 (generic function with 1 method)

Now we can finally form our image prior. For this we use a heirarchical prior where the inverse correlation length is given by a Half-Normal distribution whose peak is at zero and standard deviation is 0.1/rat where recall rat is the beam size per pixel. For the variance of the random field we use another half normal prior with standard deviation 0.1. The reason we use the half-normal priors is to prefer "simple" structures. Gaussian Markov random fields are extremly flexible models, and to prevent overfitting it is common to use priors that penalize complexity. Therefore, we want to use priors that enforce similarity to our mean image. If the data wants more complexity then it will drive us away from the prior.

cprior = HierarchicalPrior(fmap, InverseGamma(1.0, -log(0.1*rat)))
VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}(
priormap: #1
hyperprior: Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=1.5657031017929124)
θ: 0.6386906935643696
)

)

We can now form our model parameter priors. Like our other imaging examples, we use a Dirichlet prior for our image pixels. For the log gain amplitudes, we use the CalPrior which automatically constructs the prior for the given jones cache gcache.

prior = NamedDist(
         c = cprior,
         fg = Uniform(0.0, 1.0),
         σimg = truncated(Normal(0.0, 0.5); lower = 0.0),
         lgamp = CalPrior(distamp, gcache),
         gphase = CalPrior(distphase, gcachep),
        )
VLBIImagePriors.NamedDist{(:c, :fg, :σimg, :lgamp, :gphase), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Uniform{Float64}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}(
dists: (VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}(
priormap: #1
hyperprior: Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=1.5657031017929124)
θ: 0.6386906935643696
)

)
, Distributions.Uniform{Float64}(a=0.0, b=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.5); lower=0.0), CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}(
dists: DiagNormal(
dim: 129
μ: [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0]
Σ: [0.010000000000000002 0.0 … 0.0 0.0; 0.0 0.010000000000000002 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 0.010000000000000002]
)

jcache: JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}(sparse([1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [1, 1, 1, 3, 4, 4, 5, 5, 5, 7, 8, 8, 9, 9, 9, 11, 12, 12, 13, 13, 13, 15, 16, 16, 17, 17, 17, 19, 20, 20, 22, 23, 23, 24, 24, 24, 26, 27, 27, 28, 28, 28, 30, 31, 31, 32, 32, 32, 34, 35, 35, 36, 36, 36, 36, 38, 39, 39, 40, 40, 40, 41, 41, 41, 41, 43, 44, 44, 45, 45, 45, 46, 46, 46, 46, 48, 49, 49, 50, 50, 50, 51, 51, 51, 51, 51, 52, 53, 53, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 57, 57, 57, 58, 58, 59, 59, 59, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 64, 64, 64, 64, 64, 64, 65, 65, 66, 66, 66, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 71, 71, 71, 71, 71, 71, 72, 72, 73, 73, 73, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 78, 78, 78, 78, 78, 79, 79, 80, 80, 80, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 86, 86, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 95, 95, 96, 96, 98, 98, 98, 99, 100, 100, 100, 100, 100, 101, 101, 102, 102, 102, 104, 104, 104, 104, 105, 106, 106, 106, 106, 106, 107, 107, 108, 108, 108, 110, 110, 110, 110, 111, 112, 112, 112, 112, 112, 113, 113, 114, 114, 114, 116, 116, 116, 116, 117, 118, 118, 118, 118, 118, 119, 119, 120, 120, 120, 122, 122, 122, 122, 123, 124, 124, 124, 124, 124, 125, 125, 126, 126, 126, 128, 128, 128, 128, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [2, 2, 2, 3, 3, 4, 6, 6, 6, 7, 7, 8, 10, 10, 10, 11, 11, 12, 14, 14, 14, 15, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 22, 25, 25, 25, 26, 26, 27, 29, 29, 29, 30, 30, 31, 33, 33, 33, 34, 34, 35, 37, 37, 37, 37, 38, 38, 38, 39, 39, 40, 42, 42, 42, 42, 43, 43, 43, 44, 44, 45, 47, 47, 47, 47, 48, 48, 48, 49, 49, 50, 52, 52, 52, 52, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 56, 58, 58, 58, 58, 59, 59, 59, 60, 60, 60, 60, 60, 60, 61, 61, 62, 63, 63, 63, 63, 63, 65, 65, 65, 65, 66, 66, 66, 67, 67, 67, 67, 67, 67, 68, 68, 69, 70, 70, 70, 70, 70, 72, 72, 72, 72, 73, 73, 73, 74, 74, 74, 74, 74, 74, 75, 75, 76, 77, 77, 77, 77, 77, 79, 79, 79, 80, 80, 81, 81, 81, 81, 81, 82, 83, 83, 83, 83, 85, 85, 85, 86, 86, 87, 87, 87, 87, 88, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 96, 96, 97, 97, 97, 97, 98, 99, 99, 99, 101, 101, 101, 102, 102, 103, 103, 103, 103, 103, 104, 105, 105, 105, 105, 107, 107, 107, 108, 108, 109, 109, 109, 109, 109, 110, 111, 111, 111, 111, 113, 113, 113, 114, 114, 115, 115, 115, 115, 115, 116, 117, 117, 117, 117, 119, 119, 119, 120, 120, 121, 121, 121, 121, 121, 122, 123, 123, 123, 123, 125, 125, 125, 126, 126, 127, 127, 127, 127, 127, 128, 129, 129, 129, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AA, :AP, :LM, :PV, :AA, :AP, :LM, :PV, :AA, :AP  …  :AZ, :JC, :LM, :SM, :AA, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AA), (0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AA), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AA), (1.516666665673256, :AP)  …  (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AA), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), Fill(NoReference(), 25))
)
, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}(
dists: VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}(
μ: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
κ: [0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778  …  0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778]
lnorm: -180.86855628209554
)

jcache: JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}(Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([4, 5, 6, 10, 11, 12, 16, 17, 18, 22  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [2, 3, 3, 5, 6, 6, 8, 9, 9, 11, 12, 12, 14, 15, 15, 16, 17, 17, 19, 20, 20, 22, 23, 23, 25, 26, 26, 28, 29, 29, 30, 30, 30, 32, 33, 33, 34, 34, 34, 36, 37, 37, 38, 38, 38, 39, 40, 40, 42, 42, 42, 43, 43, 43, 43, 44, 44, 45, 45, 45, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 50, 50, 51, 51, 51, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 56, 56, 57, 57, 57, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 62, 62, 63, 63, 63, 65, 65, 65, 65, 66, 67, 68, 68, 70, 70, 70, 71, 71, 72, 72, 72, 74, 74, 74, 74, 75, 76, 76, 78, 78, 78, 79, 80, 80, 81, 81, 81, 83, 83, 83, 83, 84, 85, 85, 86, 86, 86, 88, 88, 88, 88, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 96, 96, 96, 98, 98, 98, 98, 99, 100, 100, 101, 101, 101, 103, 103, 103, 103, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [1, 1, 1, 2, 2, 3, 4, 4, 4, 5, 5, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 11, 11, 12, 13, 13, 13, 14, 14, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 21, 22, 22, 23, 24, 24, 24, 25, 25, 26, 27, 27, 27, 27, 28, 28, 28, 29, 29, 30, 31, 31, 31, 31, 32, 32, 32, 33, 33, 34, 35, 35, 35, 35, 36, 36, 36, 37, 37, 38, 39, 39, 39, 39, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 43, 44, 44, 44, 44, 45, 45, 45, 46, 46, 46, 46, 46, 46, 47, 47, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 52, 52, 52, 52, 52, 52, 53, 53, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 58, 58, 58, 58, 58, 58, 59, 59, 60, 61, 61, 61, 61, 61, 62, 62, 62, 63, 63, 64, 64, 64, 64, 64, 65, 66, 66, 66, 66, 67, 67, 67, 68, 68, 69, 69, 69, 69, 70, 71, 71, 71, 72, 72, 73, 73, 73, 73, 73, 74, 75, 75, 75, 75, 76, 76, 77, 77, 77, 77, 78, 79, 79, 79, 80, 80, 80, 81, 81, 82, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 85, 85, 86, 86, 87, 87, 87, 87, 87, 88, 89, 89, 89, 89, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 95, 95, 95, 96, 96, 97, 97, 97, 97, 97, 98, 99, 99, 99, 99, 100, 100, 100, 101, 101, 102, 102, 102, 102, 102, 103, 104, 104, 104, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AP, :LM, :PV, :AP, :LM, :PV, :AP, :LM, :PV, :AP  …  :AP, :AZ, :JC, :LM, :SM, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256, 1.516666665673256, 1.816666603088379  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AP), (1.516666665673256, :LM), (1.516666665673256, :PV), (1.816666603088379, :AP)  …  (7.7166666984558105, :AP), (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), SingleReference{FixedSeg{ComplexF64}}[SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AP, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))  …  SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))])
)
)
)

Putting it all together we form our likelihood and posterior objects for optimization and sampling.

lklhd = RadioLikelihood(sky, instrument, dvis; skymeta=metadata, instrumentmeta=metadata)
post = Posterior(lklhd, prior)
Posterior{RadioLikelihood{Comrade.ModelMetadata{typeof(Main.sky), NamedTuple{(:ftot, :K, :meanpr, :grid, :cache, :gcache, :gcachep), Tuple{Float64, VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}, Matrix{Float64}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}, VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}, Comrade.ModelMetadata{typeof(Main.instrument), NamedTuple{(:ftot, :K, :meanpr, :grid, :cache, :gcache, :gcachep), Tuple{Float64, VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}, Matrix{Float64}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}, VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}, Tuple{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}}, Tuple{Comrade.ConditionedLikelihood{Comrade.var"#26#27"{Float64, Vector{Float64}}, Vector{ComplexF64}}}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, NamedTuple{(:U, :V, :T, :F), NTuple{4, Vector{Float64}}}}, VLBIImagePriors.NamedDist{(:c, :fg, :σimg, :lgamp, :gphase), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Uniform{Float64}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}}(RadioLikelihood
	Number of data products: 1
, VLBIImagePriors.NamedDist{(:c, :fg, :σimg, :lgamp, :gphase), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Uniform{Float64}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}}}(
dists: (VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}(
priormap: #1
hyperprior: Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=1.5657031017929124)
θ: 0.6386906935643696
)

)
, Distributions.Uniform{Float64}(a=0.0, b=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.5); lower=0.0), CalPrior{Distributions.DiagNormal, JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}}(
dists: DiagNormal(
dim: 129
μ: [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0]
Σ: [0.010000000000000002 0.0 … 0.0 0.0; 0.0 0.010000000000000002 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 0.010000000000000002]
)

jcache: JonesCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, FillArrays.Fill{NoReference, 1, Tuple{Base.OneTo{Int64}}}}(sparse([1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [1, 1, 1, 3, 4, 4, 5, 5, 5, 7, 8, 8, 9, 9, 9, 11, 12, 12, 13, 13, 13, 15, 16, 16, 17, 17, 17, 19, 20, 20, 22, 23, 23, 24, 24, 24, 26, 27, 27, 28, 28, 28, 30, 31, 31, 32, 32, 32, 34, 35, 35, 36, 36, 36, 36, 38, 39, 39, 40, 40, 40, 41, 41, 41, 41, 43, 44, 44, 45, 45, 45, 46, 46, 46, 46, 48, 49, 49, 50, 50, 50, 51, 51, 51, 51, 51, 52, 53, 53, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 57, 57, 57, 58, 58, 59, 59, 59, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 64, 64, 64, 64, 64, 64, 65, 65, 66, 66, 66, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 71, 71, 71, 71, 71, 71, 72, 72, 73, 73, 73, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 78, 78, 78, 78, 78, 79, 79, 80, 80, 80, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 86, 86, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 95, 95, 96, 96, 98, 98, 98, 99, 100, 100, 100, 100, 100, 101, 101, 102, 102, 102, 104, 104, 104, 104, 105, 106, 106, 106, 106, 106, 107, 107, 108, 108, 108, 110, 110, 110, 110, 111, 112, 112, 112, 112, 112, 113, 113, 114, 114, 114, 116, 116, 116, 116, 117, 118, 118, 118, 118, 118, 119, 119, 120, 120, 120, 122, 122, 122, 122, 123, 124, 124, 124, 124, 124, 125, 125, 126, 126, 126, 128, 128, 128, 128, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [2, 2, 2, 3, 3, 4, 6, 6, 6, 7, 7, 8, 10, 10, 10, 11, 11, 12, 14, 14, 14, 15, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 22, 25, 25, 25, 26, 26, 27, 29, 29, 29, 30, 30, 31, 33, 33, 33, 34, 34, 35, 37, 37, 37, 37, 38, 38, 38, 39, 39, 40, 42, 42, 42, 42, 43, 43, 43, 44, 44, 45, 47, 47, 47, 47, 48, 48, 48, 49, 49, 50, 52, 52, 52, 52, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 56, 58, 58, 58, 58, 59, 59, 59, 60, 60, 60, 60, 60, 60, 61, 61, 62, 63, 63, 63, 63, 63, 65, 65, 65, 65, 66, 66, 66, 67, 67, 67, 67, 67, 67, 68, 68, 69, 70, 70, 70, 70, 70, 72, 72, 72, 72, 73, 73, 73, 74, 74, 74, 74, 74, 74, 75, 75, 76, 77, 77, 77, 77, 77, 79, 79, 79, 80, 80, 81, 81, 81, 81, 81, 82, 83, 83, 83, 83, 85, 85, 85, 86, 86, 87, 87, 87, 87, 88, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 96, 96, 97, 97, 97, 97, 98, 99, 99, 99, 101, 101, 101, 102, 102, 103, 103, 103, 103, 103, 104, 105, 105, 105, 105, 107, 107, 107, 108, 108, 109, 109, 109, 109, 109, 110, 111, 111, 111, 111, 113, 113, 113, 114, 114, 115, 115, 115, 115, 115, 116, 117, 117, 117, 117, 119, 119, 119, 120, 120, 121, 121, 121, 121, 121, 122, 123, 123, 123, 123, 125, 125, 125, 126, 126, 127, 127, 127, 127, 127, 128, 129, 129, 129, 129], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 129), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AA, :AP, :LM, :PV, :AA, :AP, :LM, :PV, :AA, :AP  …  :AZ, :JC, :LM, :SM, :AA, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AA), (0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AA), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AA), (1.516666665673256, :AP)  …  (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AA), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), Fill(NoReference(), 25))
)
, CalPrior{VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}, JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}}(
dists: VLBIImagePriors.DiagonalVonMises{Vector{Float64}, Vector{Float64}, Float64}(
μ: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
κ: [0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778  …  0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778, 0.10132118364233778]
lnorm: -180.86855628209554
)

jcache: JonesCache{Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}, NamedTuple{(:AA, :AP, :AZ, :JC, :LM, :PV, :SM), NTuple{7, ScanSeg{false}}}, Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}, Vector{SingleReference{FixedSeg{ComplexF64}}}}(Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([4, 5, 6, 10, 11, 12, 16, 17, 18, 22  …  275, 276, 277, 278, 279, 280, 281, 282, 283, 284], [2, 3, 3, 5, 6, 6, 8, 9, 9, 11, 12, 12, 14, 15, 15, 16, 17, 17, 19, 20, 20, 22, 23, 23, 25, 26, 26, 28, 29, 29, 30, 30, 30, 32, 33, 33, 34, 34, 34, 36, 37, 37, 38, 38, 38, 39, 40, 40, 42, 42, 42, 43, 43, 43, 43, 44, 44, 45, 45, 45, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 50, 50, 51, 51, 51, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 56, 56, 57, 57, 57, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 62, 62, 63, 63, 63, 65, 65, 65, 65, 66, 67, 68, 68, 70, 70, 70, 71, 71, 72, 72, 72, 74, 74, 74, 74, 75, 76, 76, 78, 78, 78, 79, 80, 80, 81, 81, 81, 83, 83, 83, 83, 84, 85, 85, 86, 86, 86, 88, 88, 88, 88, 89, 90, 90, 91, 91, 91, 93, 93, 93, 93, 94, 95, 95, 96, 96, 96, 98, 98, 98, 98, 99, 100, 100, 101, 101, 101, 103, 103, 103, 103, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), Comrade.AffineDesignMatrix{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{ComplexF64}}(sparse([2, 4, 6, 1, 5, 3, 7, 10, 12, 9  …  274, 276, 279, 283, 284, 272, 271, 275, 278, 281], [1, 1, 1, 2, 2, 3, 4, 4, 4, 5, 5, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 11, 11, 12, 13, 13, 13, 14, 14, 15, 16, 18, 18, 18, 19, 19, 20, 21, 21, 21, 22, 22, 23, 24, 24, 24, 25, 25, 26, 27, 27, 27, 27, 28, 28, 28, 29, 29, 30, 31, 31, 31, 31, 32, 32, 32, 33, 33, 34, 35, 35, 35, 35, 36, 36, 36, 37, 37, 38, 39, 39, 39, 39, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 43, 44, 44, 44, 44, 45, 45, 45, 46, 46, 46, 46, 46, 46, 47, 47, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 52, 52, 52, 52, 52, 52, 53, 53, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 57, 57, 57, 58, 58, 58, 58, 58, 58, 59, 59, 60, 61, 61, 61, 61, 61, 62, 62, 62, 63, 63, 64, 64, 64, 64, 64, 65, 66, 66, 66, 66, 67, 67, 67, 68, 68, 69, 69, 69, 69, 70, 71, 71, 71, 72, 72, 73, 73, 73, 73, 73, 74, 75, 75, 75, 75, 76, 76, 77, 77, 77, 77, 78, 79, 79, 79, 80, 80, 80, 81, 81, 82, 82, 82, 82, 82, 83, 84, 84, 84, 84, 85, 85, 85, 86, 86, 87, 87, 87, 87, 87, 88, 89, 89, 89, 89, 90, 90, 90, 91, 91, 92, 92, 92, 92, 92, 93, 94, 94, 94, 94, 95, 95, 95, 96, 96, 97, 97, 97, 97, 97, 98, 99, 99, 99, 99, 100, 100, 100, 101, 101, 102, 102, 102, 102, 102, 103, 104, 104, 104, 104], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 284, 104), ComplexF64[0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im  …  0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]), (AA = ScanSeg{false}(), AP = ScanSeg{false}(), AZ = ScanSeg{false}(), JC = ScanSeg{false}(), LM = ScanSeg{false}(), PV = ScanSeg{false}(), SM = ScanSeg{false}()), Comrade.GainSchema{Vector{Symbol}, Vector{Float64}, Vector{Tuple{Float64, Symbol}}}([:AP, :LM, :PV, :AP, :LM, :PV, :AP, :LM, :PV, :AP  …  :AP, :AZ, :JC, :LM, :SM, :AP, :AZ, :JC, :LM, :SM], [0.9166666567325592, 0.9166666567325592, 0.9166666567325592, 1.2166666388511658, 1.2166666388511658, 1.2166666388511658, 1.516666665673256, 1.516666665673256, 1.516666665673256, 1.816666603088379  …  7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.7166666984558105, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905, 7.983333349227905], [(0.9166666567325592, :AP), (0.9166666567325592, :LM), (0.9166666567325592, :PV), (1.2166666388511658, :AP), (1.2166666388511658, :LM), (1.2166666388511658, :PV), (1.516666665673256, :AP), (1.516666665673256, :LM), (1.516666665673256, :PV), (1.816666603088379, :AP)  …  (7.7166666984558105, :AP), (7.7166666984558105, :AZ), (7.7166666984558105, :JC), (7.7166666984558105, :LM), (7.7166666984558105, :SM), (7.983333349227905, :AP), (7.983333349227905, :AZ), (7.983333349227905, :JC), (7.983333349227905, :LM), (7.983333349227905, :SM)]), SingleReference{FixedSeg{ComplexF64}}[SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AP, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))  …  SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im)), SingleReference{FixedSeg{ComplexF64}}(:AA, FixedSeg{ComplexF64}(1.0 + 0.0im))])
)
)
)
)

Reconstructing the Image and Instrument Effects

To sample from this posterior, it is convenient to move from our constrained parameter space to an unconstrained one (i.e., the support of the transformed posterior is (-∞, ∞)). This is done using the asflat function.

tpost = asflat(post)
ndim = dimension(tpost)
1364

Our Posterior and TransformedPosterior objects satisfy the LogDensityProblems interface. This allows us to easily switch between different AD backends and many of Julia's statistical inference packages use this interface as well.

using LogDensityProblemsAD
using Zygote
gtpost = ADgradient(Val(:Zygote), tpost)
x0 = randn(rng, ndim)
LogDensityProblemsAD.logdensity_and_gradient(gtpost, x0)
(-1.0887331545999157e8, [-2.158201503029121, -9.600833866717904, 8.867232692624132, 0.4622198854198183, 6.233313456569947, 288.1685643097863, 6.7358488956835245, 7.67301955405554, 3.4295544290806586, 4.067083747983911  …  3666.7681059308334, 11816.75911472515, 1587.1490272240208, -343.3074914007036, 6760.359870612455, 349.8400750594293, 17534.646730431534, -7352.639076547881, -181.3718982383253, 410.6393309543062])

We can now also find the dimension of our posterior or the number of parameters we are going to sample.

Warning

This can often be different from what you would expect. This is especially true when using angular variables where we often artificially increase the dimension of the parameter space to make sampling easier.

To initialize our sampler we will use optimize using LBFGS

using ComradeOptimization
using OptimizationOptimJL
f = OptimizationFunction(tpost, Optimization.AutoZygote())
prob = Optimization.OptimizationProblem(f, prior_sample(rng, tpost), nothing)
ℓ = logdensityof(tpost)
sol = solve(prob, LBFGS(), maxiters=4_000, g_tol=1e-1);

Now transform back to parameter space

xopt = transform(tpost, sol.u)
(c = (params = [0.0010961153490788864 0.0017263185334439957 … -0.00014561030192408335 0.00019249787511360864; 0.00043755496269210685 0.0010308409108939207 … 5.029516928499581e-5 -0.00017409000762303554; … ; 0.00019548225653753936 0.000408041559832034 … -0.0002438184984174647 -0.00032739165913701543; 0.0009636703048394227 0.0013011104811493331 … -0.00015872501292635424 0.0001591992444471786], hyperparams = 0.05323026315353408), fg = 0.41770106408584806, σimg = 3.581649310814465, lgamp = [-0.0161483656442454, 0.016070087783461362, -0.8496114939704238, -0.005480487402216438, -0.0017545745019954698, 0.002037553027446116, -0.2804161575535802, 0.10112119224366847, -0.01968269074196344, 0.019755958715459068  …  -0.04792155352649273, -0.0006653893418227563, -0.6413756905942436, -0.004664363068285671, -0.005637240983404715, 0.006027328432151213, -0.05966548830727049, 0.0007476129555915554, -0.6968956900544706, -0.003228559528562449], gphase = [-0.814273784204204, 1.3889466208144312, 1.0737799998942372, -0.8711725698614121, 1.4471968698372923, 0.8732731789542477, -0.9285728943054743, 1.4867852316499277, 0.6773604953708107, -0.9773735641111784  …  -0.5846401858692766, -2.613245357131321, -2.727672355718321, -2.5968960778280126, -0.9562295805713551, -0.523845054875055, -2.7196690520749414, -2.614408159895729, -2.7602136191999316, -0.82517903903886])
Warning

Fitting gains tends to be very difficult, meaning that optimization can take a lot longer. The upside is that we usually get nicer images.

First we will evaluate our fit by plotting the residuals

using Plots
residual(vlbimodel(post, xopt), dvis)
Example block output

These look reasonable, although there may be some minor overfitting. This could be improved in a few ways, but that is beyond the goal of this quick tutorial. Plotting the image, we see that we have a much cleaner version of the closure-only image from Imaging a Black Hole using only Closure Quantities.

import WGLMakie as CM
img = intensitymap(skymodel(post, xopt), fovx, fovy, 128, 128)
imageviz(img)

Because we also fit the instrument model, we can inspect their parameters. To do this, Comrade provides a caltable function that converts the flattened gain parameters to a tabular format based on the time and its segmentation.

gt = Comrade.caltable(gcachep, xopt.gphase)
plot(gt, layout=(3,3), size=(600,500))
Example block output

The gain phases are pretty random, although much of this is due to us picking a random reference station for each scan.

Moving onto the gain amplitudes, we see that most of the gain variation is within 10% as expected except LMT, which has massive variations.

gt = Comrade.caltable(gcache, exp.(xopt.lgamp))
plot(gt, layout=(3,3), size=(600,500))
Example block output

To sample from the posterior, we will use HMC, specifically the NUTS algorithm. For information about NUTS, see Michael Betancourt's notes.

Note

For our metric, we use a diagonal matrix due to easier tuning

However, due to the need to sample a large number of gain parameters, constructing the posterior is rather time-consuming. Therefore, for this tutorial, we will only do a quick preliminary run, and any posterior inferences should be appropriately skeptical.

using ComradeAHMC
metric = DiagEuclideanMetric(ndim)
chain, stats = sample(rng, post, AHMC(;metric, autodiff=Val(:Zygote)), 700; nadapts=500, init_params=xopt)
(NamedTuple{(:c, :fg, :σimg, :lgamp, :gphase), Tuple{NamedTuple{(:params, :hyperparams), Tuple{Matrix{Float64}, Float64}}, Float64, Float64, Vector{Float64}, Vector{Float64}}}[(c = (params = [-0.00047964322683784145 0.009205483938507635 … 0.01160292510138299 -0.006834591845100404; 0.007635846551151073 0.014461713269747784 … 0.00022023919144125201 -0.003484172225582612; … ; -0.005476439080899154 0.004858228677648686 … -0.00479910350319373 0.0013000235121199187; -0.008754734102455594 0.01186583619071594 … -0.003725200156445666 -0.003159906894812609], hyperparams = 0.05329999649480524), fg = 0.41966465142142056, σimg = 3.595723578179815, lgamp = [-0.015470412234844813, 0.01576991430599477, -0.8534652120536451, -0.005688614720342752, 0.006978645611379842, -0.005893852446412423, -0.28877916261728587, 0.09182459258295772, -0.018814064850541285, 0.020157750476299828  …  -0.04935003274651614, 0.0018249223767025276, -0.6475349946907756, -0.00017967223935420835, -0.007270000468674841, 0.008006827765972951, -0.06182901595331728, -0.006620907376435598, -0.6998086600268694, 0.004493533279599257], gphase = [-0.8173370768093514, 1.3722057558793166, 1.0784708025020437, -0.868626633931342, 1.4455519819722678, 0.8708619381863993, -0.9280897735968049, 1.4870806972428716, 0.6750225467459224, -0.9787538677639246  …  -0.5869908775980214, -2.6120614352031213, -2.7292285410217008, -2.6007619918508675, -0.9615539163653786, -0.5218519918117878, -2.72018314036493, -2.6257527532001355, -2.7586500051496143, -0.8346829930951537]), (c = (params = [-0.00047964322683784145 0.009205483938507635 … 0.01160292510138299 -0.006834591845100404; 0.007635846551151073 0.014461713269747784 … 0.00022023919144125201 -0.003484172225582612; … ; -0.005476439080899154 0.004858228677648686 … -0.00479910350319373 0.0013000235121199187; -0.008754734102455594 0.01186583619071594 … -0.003725200156445666 -0.003159906894812609], hyperparams = 0.05329999649480524), fg = 0.41966465142142056, σimg = 3.595723578179815, lgamp = [-0.015470412234844813, 0.01576991430599477, -0.8534652120536451, -0.005688614720342752, 0.006978645611379842, -0.005893852446412423, -0.28877916261728587, 0.09182459258295772, -0.018814064850541285, 0.020157750476299828  …  -0.04935003274651614, 0.0018249223767025276, -0.6475349946907756, -0.00017967223935420835, -0.007270000468674841, 0.008006827765972951, -0.06182901595331728, -0.006620907376435598, -0.6998086600268694, 0.004493533279599257], gphase = [-0.8173370768093514, 1.3722057558793166, 1.0784708025020437, -0.868626633931342, 1.4455519819722678, 0.8708619381863993, -0.9280897735968049, 1.4870806972428716, 0.6750225467459224, -0.9787538677639246  …  -0.5869908775980214, -2.6120614352031213, -2.7292285410217008, -2.6007619918508675, -0.9615539163653786, -0.5218519918117878, -2.72018314036493, -2.6257527532001355, -2.7586500051496143, -0.8346829930951537]), (c = (params = [0.038741232122462825 -0.03501904497913044 … -0.11244759676503413 -0.06882054077717399; 0.1828055456702078 -0.12056856450200991 … 0.30184486674019034 0.19170207512541138; … ; -0.28407460539500645 -0.07050325259668441 … -0.1528214673116185 0.1488055218871421; -0.03494949596551905 0.3081837384993815 … -0.33201207042879877 -0.22765568372030567], hyperparams = 0.03993723586116079), fg = 0.4209387272793109, σimg = 3.671739594297515, lgamp = [-0.02236674829544234, 0.02387764719859222, -0.836321597103909, 0.023832247437366758, 0.01846558518500214, -0.019626494659254125, -0.3047704466188476, 0.09975933362834916, -0.02031356401469201, 0.020382275139787975  …  -0.02411099955314395, -0.007921159652478568, -0.6237478030052076, 0.00840219490102892, -0.01124055778086235, 0.016093311728856016, -0.06463944304941656, 0.0054104482823382136, -0.6961352607046445, -0.007692050114275841], gphase = [-0.8125939186772868, 1.3894070769109566, 1.1213237022487015, -0.8753384727130447, 1.468225841325751, 0.9244085586042193, -0.928901494396054, 1.5036462440394343, 0.7204656058243415, -0.9799334192419784  …  -0.5875717687412167, -2.6040391261711093, -2.6979983939747103, -2.573593031919027, -0.9231560988512743, -0.5317075207557125, -2.7042101528202704, -2.5868356084177444, -2.7414351614800294, -0.7982777288171875]), (c = (params = [0.08656623802355794 0.4175359719171143 … 0.2909589622884898 -0.17172872912938278; 0.20195124450420687 -0.5233782044625285 … 0.5126949217264258 0.3895709361259298; … ; 0.11699561969700426 0.16829507759356366 … -0.16499540456982537 -0.20550605944432512; -0.26194032637982295 -0.15590563835163604 … -0.271874378432749 -0.5530386022455608], hyperparams = 0.04147551299846562), fg = 0.37402121264262767, σimg = 3.4832243412611033, lgamp = [-0.014835370611388246, 0.012889197859447013, -0.7080907573950005, -0.07824923158276038, 0.0033430920544708347, 0.0002689624538821986, -0.18973970586017694, 0.03498349331968267, -0.007330187006039387, 0.007054850090507976  …  -0.11189050465134787, 0.005678433456736702, -0.6795824943310879, -0.004898919047985538, 0.010287645783438605, -0.010023695539399635, -0.12938511345489323, -0.009221801221328967, -0.7023848313187807, -0.003712691174367201], gphase = [-0.8162759094462477, 1.4760454699353178, 1.2262421296929775, -0.8680879158371378, 1.554870602659529, 1.02597984983122, -0.928906630517918, 1.5810715451389525, 0.8094560874559087, -0.9728523753023807  …  -0.5877780317266286, -2.6028667675299366, -2.740790964154436, -2.560939916908148, -0.9665393916705671, -0.5243500304012865, -2.694681552896509, -2.6147577468360446, -2.706869229759769, -0.8346561091240836]), (c = (params = [0.17568096039809955 0.02614232819626226 … -0.10698276641752319 -0.3755104642963316; 0.2952640212580534 0.2122163758131522 … 0.7574508536451688 -0.04292855439184173; … ; 0.32883965643412305 -0.163399524193005 … -0.5241497347384707 0.15000160188830788; -0.12912049016213145 -0.5482934976008565 … 0.22388862283312808 -0.389586294042736], hyperparams = 0.16433230571950402), fg = 0.46579937801069105, σimg = 2.8337735814917844, lgamp = [-0.0220334985130115, 0.023959127793800945, -0.8495891223090004, 0.04714133307912389, -0.009277714776368323, 0.005393078370545133, -0.27355132062374216, 0.14181682373269502, -0.007738608291530142, 0.009567882613080189  …  0.01013048921433802, -0.014609271127410305, -0.5351922927617322, 0.0014390726660981476, -0.003831569857497551, 0.002117597471702082, -0.011477838518220328, -0.007464547955012769, -0.5994976097656235, -0.0050009128821801654], gphase = [-0.81844132881122, 1.4090462888689785, 1.1228317222327724, -0.8748677752883902, 1.4806293673318733, 0.9194841365996675, -0.9289047234688206, 1.5277510940878831, 0.7416349992766103, -0.9808701589291358  …  -0.5898567963806869, -2.5896864394121204, -2.9151703354161764, -2.605256958290375, -1.138765268057353, -0.5181818781136914, -2.689154250899593, -2.8037869067686447, -2.77894140256522, -1.0098211003557078]), (c = (params = [0.3556080278462512 0.058341658325282104 … -0.012067360283614229 -0.21024793702932681; 0.2086503699594264 0.3537621024738001 … 0.7342404845954449 0.12494483855466268; … ; 0.400787303068321 -0.13178733022983324 … -0.37899960949497763 0.24399214617641407; -0.008275330430147694 -0.4843123872838715 … 0.15286051027231895 -0.3670557490256242], hyperparams = 0.16324714022638165), fg = 0.4735607608516858, σimg = 2.760253206112478, lgamp = [-0.003919507566381899, 0.0018920738328059503, -0.8313852707512406, 0.009392772049916774, 0.005836829733825127, 0.0006411493709764179, -0.23330490726007852, 0.11737021086532844, -0.01333345868935729, 0.0101449951194571  …  0.019354736536856994, -0.001429216503523952, -0.5438039724368028, 0.0006552162758944302, -0.009474361720986537, 0.010670681951379307, 0.014691192761817374, -0.006922361969674288, -0.5942251049195979, -0.00716441920015399], gphase = [-0.8083464515619065, 1.4078270187207882, 1.1410465372565703, -0.867741244138999, 1.4385980650665084, 0.9434167724587523, -0.9307031687347945, 1.500039493765092, 0.740210660834176, -0.9773762242642431  …  -0.5814310731025978, -2.591430105206676, -2.9047926466899376, -2.5865657028058466, -1.1285470444000527, -0.5269076273362062, -2.6965680192860324, -2.8093399308506606, -2.749764780752682, -1.0137427017290153]), (c = (params = [0.6330007459693808 0.900872659185292 … 0.3616390939891469 -0.06703361625427014; -0.4343149944197584 0.26020841345643797 … -0.025500326880152687 0.1393191034112235; … ; -0.09603514060888312 0.4040516538608081 … -0.2224975373360066 -1.093621658411774; 0.6679489622183911 -0.6731615041080937 … -0.8657785420964272 -0.1326310002032061], hyperparams = 0.1450629930498135), fg = 0.4203454656753338, σimg = 2.587989410797412, lgamp = [-0.0035218605955253594, 0.0027812019748388655, -0.8735787709547135, -0.016068984409517744, 0.013720848056292215, -0.0007737784144117483, -0.28837030688410004, 0.10558988280069725, -0.021750518552181294, 0.019880846833009863  …  -0.030967944006604894, -0.0006917114581427673, -0.7425947479037593, -0.013800802684580503, -0.009904336801389073, 0.0009506867003550047, -0.02750054602941346, -0.012664666513020154, -0.8055789406844872, -0.006985514815560444], gphase = [-0.8171278878429584, 1.4290952920295878, 1.1545549118773593, -0.8733013444976081, 1.4390125877608388, 0.9445942005623701, -0.9303120680339587, 1.481447474693052, 0.7476571404468336, -0.9774683865279958  …  -0.5829585679579111, -2.6670359366606484, -2.9406309626130467, -2.5259101200667895, -1.1538287047074682, -0.5208809561914081, -2.759892231247331, -2.8309327118271996, -2.6961846206852704, -1.0335826789100704]), (c = (params = [0.27725210761545643 0.19764031264778478 … -0.020310369995121143 0.034507821019764164; 0.29229126632596186 -0.05313500780724969 … 0.6063398083081618 0.7397882318187567; … ; -0.6686417977470489 -0.2141439322646931 … -1.109653054031386 -0.8814344995555915; -0.5370316965793925 0.2851621018930625 … -0.39700756821215316 -0.7993229267289136], hyperparams = 0.21467818559839896), fg = 0.4389298309131554, σimg = 2.719616192888128, lgamp = [-0.04163623144221909, 0.04025491002064998, -0.8245862995178127, -0.002441673395561491, -0.0028454037515815373, 0.0073581114615081, -0.2723887616880537, 0.09509128867767072, -0.015957355668338565, 0.019929770904518248  …  -0.025676870268877952, 0.008151139951705773, -0.6378366016092765, 0.0032489170277324618, -0.0032005785854518996, 0.007211566125483422, -0.07163684701383476, 0.011568663255825076, -0.6997006564157452, -0.017284457991005683], gphase = [-0.8067730610748859, 1.4632408658518097, 1.0792777752347342, -0.8727479129759071, 1.4681450373378742, 0.8772351424495671, -0.9268989717372862, 1.5087924600880203, 0.6786710887665395, -0.9768613204341476  …  -0.5900274645107454, -2.7071340941883197, -2.948472644950028, -2.5645505433725493, -1.1685627402541308, -0.5249203243195507, -2.8082468068183477, -2.7691582435237185, -2.733347167147591, -0.9802521365615637]), (c = (params = [0.39220405496199007 0.14899871167004614 … 0.6039742349906062 -0.1916772327290916; 0.5686500982495136 0.03370086615811059 … 0.7949407751320374 0.6291962843873322; … ; -0.8094107692230892 -0.4731670211547474 … -0.9577546384596477 -1.618901458657606; -0.7568612148398316 0.6778537729295132 … -1.0099741083818385 -0.2821985783359817], hyperparams = 0.37174861283280025), fg = 0.4163615264479086, σimg = 2.370189147462092, lgamp = [-0.01819741799033549, 0.02262366037561621, -0.8571226319886492, 0.00817119179772234, 0.006586689080619326, -0.006636014456251399, -0.3056539453421313, 0.10823530808824973, 0.005369294834140274, -0.004774682631058243  …  -0.04436969484249259, -0.011191081632016958, -0.6667900413593256, -0.02560089135432416, -0.019682665510362338, 0.0213109072495775, -0.04924650274688373, 0.007096754668046719, -0.7372342220805242, -0.011513243198868043], gphase = [-0.8152545642200215, 1.4027870250603016, 1.2084834380563332, -0.872547543326057, 1.4219925488704004, 0.969144060825975, -0.9263217830396068, 1.4614660662464845, 0.7433197510314935, -0.9747965729074944  …  -0.5797235042061338, -2.7320986519957877, -2.9900131605597275, -2.597373683713413, -1.2299897494207932, -0.524618737425895, -2.845004492827115, -2.8518743649200706, -2.7762201057339464, -1.0626201611102022]), (c = (params = [0.6126984779114261 0.3626972099266199 … 0.42547111633900053 -0.3108952659437649; 0.4927187565240049 0.10690750399420203 … 0.52471226193724 0.5363578182022758; … ; -0.6707640407630969 -0.621996823009847 … -0.8016216193597239 -1.5221217260674142; -0.6607152793372868 0.7595250645370911 … -0.932567143445629 -0.6241012930398498], hyperparams = 0.35814028479734733), fg = 0.3883543027447048, σimg = 2.364802047558366, lgamp = [0.011251657388250492, -0.017341534143241227, -0.9511558848312075, -0.04989765687308729, -0.018711029591788858, 0.01754051920130204, -0.3389630250253764, 0.087229265524659, -0.024891166595177094, 0.025368610566702504  …  -0.05218292064154218, 0.0024446961797397505, -0.6973745629209211, -0.022392771539672614, -0.016572674885176664, 0.01843690208340859, -0.046185125618991456, 0.007554746534235415, -0.7667997799769319, -0.005540310288765093], gphase = [-0.810598901921097, 1.4062834146900116, 1.1377055394966267, -0.8722732590475004, 1.4122078555479378, 0.924356738991154, -0.9277653381868372, 1.45095746169178, 0.7257398678954357, -0.9780394292836179  …  -0.5789324675923788, -2.7232066787171316, -3.0546097064882796, -2.6306123383210913, -1.2772878394095395, -0.5214386139488802, -2.8133470725456107, -2.8630763317137227, -2.776294665316427, -1.0776521842560158])  …  (c = (params = [1.3309256393485462 0.29188346335234 … 1.2967307775143881 1.7106358343942496; 0.4621398890422251 0.14963344763218484 … 1.3265768302056493 0.22444191636013797; … ; 1.2567961809050985 1.86993158262528 … 1.6756403728679676 1.603264319396258; 1.4842524070716285 2.0482150850115275 … 1.568874671517072 2.45201754251245], hyperparams = 3.1656949234219605), fg = 0.463894122829631, σimg = 2.3541474538522658, lgamp = [-0.01918029868818812, 0.01732263723679625, -0.8692692070967029, -0.04652007941782835, -0.002233200760971143, 0.005142371047870665, -0.33268731415885033, 0.07782944655777373, -0.03253947464731174, 0.03174221712102438  …  0.011963068154569002, -0.007907630436211578, -0.8653405950456227, 0.00011468387983423365, -0.030734404625474634, 0.029033284791606315, 0.005934312285312915, 0.004473136534318901, -0.9263558113570702, 0.0048637566862855764], gphase = [-0.8108740338541437, -3.1234893151052927, 2.466455894072509, -0.872823498192551, -3.009555786328437, 2.277866478640862, -0.9318271947966507, -2.956623456334878, 2.075107806592947, -0.984077981375479  …  -0.5892252865326775, -0.6100232693576094, -0.4564345996354273, -0.7682485429414002, 1.313950449193838, -0.525076099367594, -0.7328792829904172, -0.426538845729417, -0.9513988292461102, 1.3763135422770867]), (c = (params = [1.453277060386724 0.28131146412944213 … 1.7119573812494768 1.6871468137883194; 0.14330050065933603 0.4918339004348529 … 0.740589659904402 -0.034615315503987436; … ; 1.6790989302494586 1.906183097003636 … 1.4899275057360593 2.1475637144592015; 1.2050618602617817 1.5061783441177434 … 1.8994402909082824 1.7651750491766247], hyperparams = 2.680010963284125), fg = 0.4541229825918648, σimg = 2.424989223707113, lgamp = [-0.02615774855932059, 0.03239243089081475, -0.7562821389594119, -0.013645120507640562, 0.0437053550566146, -0.04840490407419193, -0.24566022532701257, 0.04620419160350033, -0.009292114455855887, 0.007623490799461198  …  -0.022828767317210678, -0.0007128985849700902, -0.8026485540887128, -0.01361286366973952, 0.011126353686137711, -0.008012974292259445, -0.027798866075858367, -0.0024108891647277647, -0.8498094488822678, 0.0034002219146179696], gphase = [-0.8170863478963206, -3.029271475690606, 2.5055749391111597, -0.8737284332644093, -2.9683783241627237, 2.3100516141705314, -0.9273989996174119, -2.9244702866864727, 2.112915561598558, -0.9753905064502324  …  -0.5898959536459051, -0.6317983288975199, -0.3756845353852381, -0.7494132205254576, 1.3862455361393784, -0.5241701834798975, -0.7562759653356387, -0.32239686249841565, -0.9544238304567628, 1.4816559457325065]), (c = (params = [1.1653008526296202 0.3716373876269395 … 1.7714362116681253 1.7572185619340612; 0.04310753796430683 0.5405957993269281 … 0.720268523538502 0.1376697551858162; … ; 1.4153394423242334 1.7853036226356815 … 1.6039629142951075 2.147512641635318; 1.3049003582226146 1.5912705007141528 … 1.8760682238140463 1.8529971722176966], hyperparams = 2.334297276461947), fg = 0.4616884253289024, σimg = 2.5081413813296374, lgamp = [0.009041658899414344, -0.014744923278313431, -0.7518534225601115, -0.027635901047464603, 0.03634687910083518, -0.03858744594878112, -0.1950630480647946, 0.0827256049480983, -0.026518398820742795, 0.026884080682003764  …  -0.005002252805611547, -0.006298929848932885, -0.7797540405224306, -0.0018953532696786202, 0.004549419005974079, 0.0030307551659256617, -0.044328619522677554, 0.007776068253341058, -0.843530182282263, 0.003268709672495086], gphase = [-0.8112374891594419, -3.0314577569742682, 2.4850277327882484, -0.8736595940597502, -3.0085490665448087, 2.2802618366685725, -0.9271024157751504, -2.9501311385289406, 2.0924507633080838, -0.97689889121984  …  -0.5839996624359597, -0.6386483879579232, -0.38922940363153674, -0.7833545740814302, 1.3863445826282796, -0.5226026912725009, -0.7749251371593293, -0.325194109255075, -0.9900911604207122, 1.4507163156762364]), (c = (params = [1.0915390303936372 0.6052612682867758 … 1.4232008263907932 1.3196767663395053; 0.42703581531354307 0.39922547105244655 … 0.36616559339233595 0.3359785028548655; … ; 1.178734390350018 1.856883546126555 … 1.53056166734284 1.9881351252127408; 1.592683838073868 1.3360218666844699 … 1.755525541914787 1.9855826718045517], hyperparams = 3.038740250666408), fg = 0.4776631155341027, σimg = 2.5063875522508265, lgamp = [-0.0233875433406629, 0.019811647793840422, -0.8513846806297128, -0.026949497251087008, -0.005625430716213517, 0.009592021850101313, -0.3076277451433557, 0.09915745857315617, -0.003336980746301324, -0.0010009766643089449  …  0.00452284217286228, 0.007129300751276383, -0.8215674750847126, -0.013972800184883156, -0.008754372204934463, 0.007373600148070171, -0.017841566849060003, 0.018374100427447127, -0.8782447859913676, -0.006248357466213856], gphase = [-0.8208398104976609, -3.0811273591784265, 2.5588931158473813, -0.866971911237119, -3.0011476986180496, 2.362280802966423, -0.9279026342324125, -2.935916849393927, 2.1666009113709768, -0.9787825819781657  …  -0.585502725150239, -0.6156296347957985, -0.4247764976336921, -0.7382420063117373, 1.349801158214739, -0.5258540590915654, -0.7518973819424795, -0.39505030457683693, -0.9491197506114055, 1.3975168905763997]), (c = (params = [1.1106930079777422 0.36179700859511366 … 0.7038884657427229 1.5708136564143855; 1.4993084950478384 0.9759115557286182 … 1.0664178095978973 0.6324444609556895; … ; 0.5964331309541305 1.6666738214924792 … 1.0293199424350565 1.244661947252239; 1.1691437510217737 0.9265698377931929 … 2.080491884370604 1.35863748025063], hyperparams = 2.1843874471931013), fg = 0.47761765391586086, σimg = 2.097316161173144, lgamp = [-0.020044413002742192, 0.015570627695603646, -0.7884323874088603, -0.06639007028352238, -0.012410419159078949, 0.011963874424308737, -0.2138592784451916, 0.0487404072276905, 0.0038894003943495165, -0.0022035065482500548  …  -0.020352780729352503, 0.00933188134502804, -0.8010391322820328, -0.005865797551669514, -0.006323414818856738, 0.009381505189405458, -0.05231018617580104, -0.005514616913496795, -0.8661730403692472, 0.005635759052981004], gphase = [-0.819372200335615, -3.0841652177215195, 2.5074729615397993, -0.8701731337473169, -3.019867045541261, 2.3109116737077646, -0.9299820914089447, -2.9808812280676094, 2.117680568319852, -0.9746800070187251  …  -0.5830385081287222, -0.6733145920609613, -0.6129439859498051, -0.7980741883769845, 1.1399558220580908, -0.521901233870705, -0.7727789125631359, -0.6024445976221513, -0.9935092147150577, 1.185870024190689]), (c = (params = [1.8602141572297353 0.9607248755947884 … 1.108433585442198 1.582459445763952; 1.8053347189417186 0.9772170430046513 … 0.6453940071821492 0.667465512857211; … ; 0.7054770233107094 1.1218655245781541 … 1.0705693934755378 1.3574532520750509; 1.4893766766533274 0.46059309178927677 … 2.1585420755900016 2.282848261041642], hyperparams = 2.974586133979504), fg = 0.5588040045688297, σimg = 2.5568398553964817, lgamp = [-0.016807402513654777, 0.025857736514460886, -0.7513677847050884, -0.01823695075200223, 0.0060318858800895445, -0.0007424196907430147, -0.1977023973686524, 0.09476916225511156, 0.018720446226549742, -0.02025444382978698  …  0.023568859861274976, -0.0007717273543095327, -0.7208525684236512, 0.001875314582826382, 0.013826628359797916, -0.007739843613291754, 0.015543516821347678, 0.006484019033691749, -0.7632404639770997, -0.005058080327532978], gphase = [-0.807974816532237, -2.9560971403877416, 2.5528983259071265, -0.8760841750749198, -2.9371385396452436, 2.3426207111953223, -0.9278187606380714, -2.893275720201029, 2.158166350197892, -0.9759999906760961  …  -0.5851119762497791, -0.6491172418035464, -0.6069701328427073, -0.758678727758537, 1.1696647992590423, -0.5215596477232184, -0.7477927510787216, -0.5378453558179277, -0.9637915999316866, 1.2285177412664914]), (c = (params = [2.2703182073544053 1.3658197161436019 … 1.7478504293092512 1.4827246233936504; 2.0426103239693227 0.7325474728207155 … 0.43751523747430643 0.4777867566980796; … ; 1.717595335494885 1.0951773291253966 … 1.3215926264627327 2.0520167342429896; 1.8294607569525783 0.6717756141585816 … 2.0579917648653363 2.5032310974362737], hyperparams = 4.112187524189807), fg = 0.39628545402247073, σimg = 2.1845560496049172, lgamp = [-0.010990471110754562, 0.005979956416975767, -0.9603098665323744, -0.006385730634761613, 0.02476530693923549, -0.02754619290398486, -0.4242223805643814, 0.10404614850701516, 0.0005475046831620992, -0.0024993579251590906  …  -0.03642424456018684, -0.002186610751007387, -0.9652839475086331, -0.007187273324431484, -0.03755106836065967, 0.030729231920844015, -0.037676522749134335, -0.0013109364766083949, -1.0252918685953494, 0.003926660264121333], gphase = [-0.8131893238637907, -3.1059329239882283, 2.619233344364562, -0.86632774719104, -3.041541412738752, 2.402074931759023, -0.9281243729314926, -3.016185550261608, 2.1956182556104022, -0.9789524606054778  …  -0.5788465645390266, -0.6816126678844805, -0.596953831075423, -0.8199010722039113, 1.175319448452959, -0.5228309165160693, -0.8117859080877099, -0.5686048892570389, -1.0287011916151214, 1.2332087536756502]), (c = (params = [2.3569957863871096 1.3757012041785912 … 1.7135713060031927 1.4463609406957512; 2.0448563521958514 1.0216117808030307 … 0.3430769951049735 0.6040974836219135; … ; 1.692758182007633 0.9243005696770777 … 1.46337921466287 2.199721609052218; 1.6319174738698807 0.6176773265453305 … 2.2507132630581625 2.4459556765918062], hyperparams = 3.275337438590304), fg = 0.4351442471904084, σimg = 2.2768875877564896, lgamp = [-0.02368514186044632, 0.022386162300110864, -0.9640963682235609, 0.012304719596531252, 0.028656052199234314, -0.02885815189987781, -0.41286560494987146, 0.113757295781278, 0.011651197461965973, -0.0065606190124679175  …  -0.027182468064736905, -0.007153248642153967, -0.9434043814888834, -0.0015695372299573715, -0.005376392181080315, -0.000997423375441779, -0.04732198550844157, 0.0024446108877889464, -1.015347033671655, 0.005921609232399299], gphase = [-0.8125632249439046, -3.071606379800115, 2.6980152315201367, -0.8666205192291729, -3.005969945519466, 2.477375458385608, -0.9281219479379491, -2.975435775305624, 2.2656686648884743, -0.9763123929734601  …  -0.5868304684932142, -0.6924936756152285, -0.5890722291239194, -0.8118087902997759, 1.1800868932647233, -0.5296768672562202, -0.8047904438741716, -0.5882639626854911, -1.0117869523523861, 1.2164982859310034]), (c = (params = [1.8690666872597845 1.1992522469836073 … 1.730760960400022 1.102401165953885; 1.8624161574732934 1.0822418514776688 … 0.49421527931943393 0.6032946618380935; … ; 0.9589403646485817 0.38688576479908854 … 1.4427851702002568 2.253613158656764; 1.8636375264503589 0.7526946800280306 … 2.2793805311117548 2.846390448619223], hyperparams = 3.8182495678555064), fg = 0.4587864898656686, σimg = 2.159276588045723, lgamp = [-0.01879243142576986, 0.019379450468412524, -0.9021539411301128, -0.05303920410387659, -0.008876102879473595, 0.007450929028388069, -0.33934622586936436, 0.08929388422026738, 0.010750513631775302, -0.018236160327707904  …  0.006202260368300759, -0.0027548293656044385, -0.8733590749193348, 0.0017301745019983293, -0.011663223280720073, 0.016822447643000234, -0.021744016050257893, -0.004076081694696158, -0.9473008828217123, -0.003966705592095135], gphase = [-0.8181518457241498, -3.024732515245587, 2.6335628517019263, -0.865242232413232, -2.9564651695747997, 2.4207562534152167, -0.9289487772751372, -2.9254403363631654, 2.213303124761485, -0.9787926932230867  …  -0.5822585824923011, -0.6510450820850435, -0.6245670335462229, -0.7851276170856785, 1.1379632092688827, -0.5288106660986361, -0.7514132834248239, -0.6016600542976622, -0.9699428584231476, 1.1820691723719747]), (c = (params = [1.7590829250916755 1.8906863271847147 … 2.0428809915767525 1.6268340535459809; 1.3249207601119224 0.23036541677064462 … 1.0106500420018671 1.341106407554896; … ; 0.6350030656916829 0.0525833266533711 … 1.6488824109334221 2.0986046172602193; 2.012725051145101 0.23507718710672384 … 3.0184462025564107 2.81478237624103], hyperparams = 4.775656596780323), fg = 0.5108289195639003, σimg = 2.4233309087692865, lgamp = [-0.0276233152261982, 0.03231362246050902, -0.8549899249020164, -0.010336290584973377, -0.004321674999294648, -0.00021297200326843898, -0.35617661892610786, 0.11807644634585489, -0.03636244411160959, 0.034636930824724244  …  0.010091559060628941, 0.0005168239749250276, -0.8397843429383898, -0.010218834802605084, -0.00997199417779513, 0.009490495593028645, 0.003302743531536194, 0.009377064998083225, -0.9119559647793445, 0.00787686556601841], gphase = [-0.8179364026739211, -2.9777884876680973, 2.668093966445167, -0.8675300626852824, -2.966488967805061, 2.4652765028316725, -0.9308894844353438, -2.9420212689793424, 2.251475716346789, -0.9752187042587872  …  -0.5879327145334473, -0.6285563042498951, -0.669803557068671, -0.7942869960244423, 1.0835714843413693, -0.5169155647121144, -0.7445832007293773, -0.6513630205259877, -0.9870207610914347, 1.1296694720069127])], NamedTuple{(:n_steps, :is_accept, :acceptance_rate, :log_density, :hamiltonian_energy, :hamiltonian_energy_error, :max_hamiltonian_energy_error, :tree_depth, :numerical_error, :step_size, :nom_step_size, :is_adapt), Tuple{Int64, Bool, Float64, Float64, Float64, Float64, Float64, Int64, Bool, Float64, Float64, Bool}}[(n_steps = 1023, is_accept = 1, acceptance_rate = 0.992088657647762, log_density = 2023.1203218859935, hamiltonian_energy = -1455.5059965742194, hamiltonian_energy_error = 0.008068993188317108, max_hamiltonian_energy_error = 0.01204978010468949, tree_depth = 10, numerical_error = 0, step_size = 0.0001, nom_step_size = 0.0001, is_adapt = 1), (n_steps = 3, is_accept = 1, acceptance_rate = 0.003418539023472627, log_density = 2023.1203218859935, hamiltonian_energy = -1369.1886620942337, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 1164.313684101026, tree_depth = 1, numerical_error = 1, step_size = 0.0015529575284482515, nom_step_size = 0.0015529575284482515, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9715069332864616, log_density = 1889.6828322465426, hamiltonian_energy = -1329.6210639563033, hamiltonian_energy_error = 0.03731784363526458, max_hamiltonian_energy_error = 0.07777409170512328, tree_depth = 10, numerical_error = 0, step_size = 0.0003044948999276269, nom_step_size = 0.0003044948999276269, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9921746544422764, log_density = 1718.2394158509326, hamiltonian_energy = -1215.701744947508, hamiltonian_energy_error = -0.04523408007526086, max_hamiltonian_energy_error = -0.14469889472229625, tree_depth = 10, numerical_error = 0, step_size = 0.0004704476349862903, nom_step_size = 0.0004704476349862903, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.5026715788982943, log_density = 1516.8765434820955, hamiltonian_energy = -1083.4797492956864, hamiltonian_energy_error = 0.9948932784595854, max_hamiltonian_energy_error = 1.771680767240241, tree_depth = 10, numerical_error = 0, step_size = 0.0008899380075297112, nom_step_size = 0.0008899380075297112, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9972920315528673, log_density = 1507.7506235021976, hamiltonian_energy = -861.2743652054917, hamiltonian_energy_error = -0.06907184293186219, max_hamiltonian_energy_error = -0.24223208973262444, tree_depth = 10, numerical_error = 0, step_size = 0.00042355732262230596, nom_step_size = 0.00042355732262230596, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.6915504266237312, log_density = 1370.9184806229496, hamiltonian_energy = -779.3307154266065, hamiltonian_energy_error = 0.6231647744625661, max_hamiltonian_energy_error = 1.262501566984156, tree_depth = 10, numerical_error = 0, step_size = 0.0008824353617909085, nom_step_size = 0.0008824353617909085, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9575140228730465, log_density = 1398.1453997002843, hamiltonian_energy = -696.0287250346626, hamiltonian_energy_error = -0.5359941842700664, max_hamiltonian_energy_error = -0.811221426633665, tree_depth = 10, numerical_error = 0, step_size = 0.0007341240759536836, nom_step_size = 0.0007341240759536836, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.24751735461939145, log_density = 1383.432525661521, hamiltonian_energy = -758.5565090375837, hamiltonian_energy_error = -0.7380830070115962, max_hamiltonian_energy_error = 849.9051247731503, tree_depth = 10, numerical_error = 0, step_size = 0.0014050778788120004, nom_step_size = 0.0014050778788120004, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9852785294889369, log_density = 1359.551206269159, hamiltonian_energy = -702.7887076227364, hamiltonian_energy_error = 0.008619578556590568, max_hamiltonian_energy_error = 0.0582255461265504, tree_depth = 10, numerical_error = 0, step_size = 0.0002879265408143687, nom_step_size = 0.0002879265408143687, is_adapt = 1)  …  (n_steps = 1023, is_accept = 1, acceptance_rate = 0.8798365951878462, log_density = 1278.602278231098, hamiltonian_energy = -607.4779697516094, hamiltonian_energy_error = -0.44809241566770197, max_hamiltonian_energy_error = -0.6768003826231279, tree_depth = 10, numerical_error = 0, step_size = 0.0008272116664030766, nom_step_size = 0.0008272116664030766, is_adapt = 0), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.7597429982626004, log_density = 1288.9310788916548, hamiltonian_energy = -590.2213649540859, hamiltonian_energy_error = 0.5113322093903889, max_hamiltonian_energy_error = 0.7027848333015072, tree_depth = 10, numerical_error = 0, step_size = 0.0008272116664030766, nom_step_size = 0.0008272116664030766, is_adapt = 0), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.8873359565425389, log_density = 1268.5834615341287, hamiltonian_energy = -620.7150225072978, hamiltonian_energy_error = 0.16344768032035972, max_hamiltonian_energy_error = 0.726605019723138, tree_depth = 10, numerical_error = 0, step_size = 0.0008272116664030766, nom_step_size = 0.0008272116664030766, is_adapt = 0), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.981498520828045, log_density = 1287.7796063390224, hamiltonian_energy = -574.5127310524122, hamiltonian_energy_error = -0.3812605560333395, max_hamiltonian_energy_error = -0.777350639879387, tree_depth = 10, numerical_error = 0, step_size = 0.0008272116664030766, nom_step_size = 0.0008272116664030766, is_adapt = 0), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.8830947738126788, log_density = 1275.999897174842, hamiltonian_energy = -591.7285065014979, hamiltonian_energy_error = 0.0680938110384659, max_hamiltonian_energy_error = 0.5294166343253437, tree_depth = 10, numerical_error = 0, step_size = 0.0008272116664030766, nom_step_size = 0.0008272116664030766, is_adapt = 0), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9878156088283871, log_density = 1262.9805964682778, hamiltonian_energy = -573.2227016081032, hamiltonian_energy_error = -0.3743285444494404, max_hamiltonian_energy_error = -0.5784581694830422, tree_depth = 10, numerical_error = 0, step_size = 0.0008272116664030766, nom_step_size = 0.0008272116664030766, is_adapt = 0), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.6933111451355685, log_density = 1290.072028289198, hamiltonian_energy = -618.409534427879, hamiltonian_energy_error = 0.41655001561321114, max_hamiltonian_energy_error = 0.8916482186190251, tree_depth = 10, numerical_error = 0, step_size = 0.0008272116664030766, nom_step_size = 0.0008272116664030766, is_adapt = 0), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9448940332720982, log_density = 1295.868489887484, hamiltonian_energy = -646.6102498585763, hamiltonian_energy_error = 0.042960902328559314, max_hamiltonian_energy_error = 0.43054163152203273, tree_depth = 10, numerical_error = 0, step_size = 0.0008272116664030766, nom_step_size = 0.0008272116664030766, is_adapt = 0), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9516494126126377, log_density = 1280.1445306447251, hamiltonian_energy = -627.3117993183696, hamiltonian_energy_error = -0.02283563012633749, max_hamiltonian_energy_error = -0.5157014278593124, tree_depth = 10, numerical_error = 0, step_size = 0.0008272116664030766, nom_step_size = 0.0008272116664030766, is_adapt = 0), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.887197226897276, log_density = 1269.6137866370314, hamiltonian_energy = -589.4713943343947, hamiltonian_energy_error = -0.4644369681500393, max_hamiltonian_energy_error = -0.7268966082738189, tree_depth = 10, numerical_error = 0, step_size = 0.0008272116664030766, nom_step_size = 0.0008272116664030766, is_adapt = 0)])
Note

The above sampler will store the samples in memory, i.e. RAM. For large models this can lead to out-of-memory issues. To fix that you can include the keyword argument saveto = DiskStore() which periodically saves the samples to disk limiting memory useage. You can load the chain using load_table(diskout) where diskout is the object returned from sample. For more information please see ComradeAHMC.

Now we prune the adaptation phase

chain = chain[501:end]
Table with 5 columns and 200 rows:
      c                     fg        σimg     lgamp                 ⋯
    ┌─────────────────────────────────────────────────────────────────
 1  │ (params = [-0.70433…  0.454517  2.15958  [0.0125397, -0.0129…  ⋯
 2  │ (params = [-0.32606…  0.501913  2.55172  [-0.0153254, 0.0167…  ⋯
 3  │ (params = [-0.80183…  0.430375  2.04608  [-0.0269995, 0.0275…  ⋯
 4  │ (params = [-0.88897…  0.460489  1.99411  [0.0144494, -0.0151…  ⋯
 5  │ (params = [-0.69288…  0.462881  2.20096  [-0.0266414, 0.0298…  ⋯
 6  │ (params = [-0.50557…  0.46705   2.45141  [-0.0050745, 0.0081…  ⋯
 7  │ (params = [-0.05399…  0.522794  2.2454   [-0.0264854, 0.0251…  ⋯
 8  │ (params = [-0.18360…  0.471547  2.41993  [0.000806183, 0.001…  ⋯
 9  │ (params = [-0.18540…  0.453225  2.3794   [-0.050117, 0.04962…  ⋯
 10 │ (params = [0.446518…  0.4547    2.16084  [0.0126505, -0.0137…  ⋯
 11 │ (params = [0.964229…  0.547282  2.371    [-0.00941108, 0.010…  ⋯
 12 │ (params = [0.913393…  0.457493  2.06352  [-0.0190426, 0.0212…  ⋯
 13 │ (params = [-0.04555…  0.540161  2.33694  [-0.0422716, 0.0379…  ⋯
 14 │ (params = [-0.10516…  0.47533   2.25502  [-0.0176551, 0.0144…  ⋯
 15 │ (params = [0.192245…  0.505185  2.41026  [-0.00140757, -0.00…  ⋯
 16 │ (params = [0.43332 …  0.478945  2.13724  [0.00191063, -0.000…  ⋯
 17 │ (params = [0.567683…  0.488336  2.24032  [0.0028064, -0.0052…  ⋯
 ⋮  │          ⋮               ⋮         ⋮              ⋮            ⋱
Warning

This should be run for likely an order of magnitude more steps to properly estimate expectations of the posterior

Now that we have our posterior, we can put error bars on all of our plots above. Let's start by finding the mean and standard deviation of the gain phases

gphase  = hcat(chain.gphase...)
mgphase = mean(gphase, dims=2)
sgphase = std(gphase, dims=2)
104×1 Matrix{Float64}:
 0.003976806014595015
 2.8139705247950726
 0.44953293577540865
 0.0032016328193606356
 2.7850961932535845
 0.456037412271608
 0.0031349781767195314
 2.765865711760177
 0.4629786795079212
 0.003651918633545869
 ⋮
 0.3404429495865237
 0.3096630073459648
 0.27459907977070996
 0.30972301735227287
 0.004255636202272551
 0.3382628525539642
 0.2984338631891679
 0.27167323197735216
 0.29754939139840036

and now the gain amplitudes

gamp  = exp.(hcat(chain.lgamp...))
mgamp = mean(gamp, dims=2)
sgamp = std(gamp, dims=2)
129×1 Matrix{Float64}:
 0.02200554378173898
 0.02237863834527218
 0.03498923610726912
 0.04249374713781529
 0.017410287315397226
 0.017537750173159342
 0.04857677043418206
 0.042575790108844255
 0.016262360550605332
 0.016844173696842046
 ⋮
 0.009813160094777248
 0.02083774219071793
 0.009577719930036007
 0.014508609293675642
 0.014825868857654976
 0.03408823551014276
 0.009491831673921307
 0.020161505446656083
 0.009369986311727285

Now we can use the measurements package to automatically plot everything with error bars. First we create a caltable the same way but making sure all of our variables have errors attached to them.

using Measurements
gmeas_am = measurement.(mgamp, sgamp)
ctable_am = caltable(gcache, vec(gmeas_am)) # caltable expects gmeas_am to be a Vector
gmeas_ph = measurement.(mgphase, sgphase)
ctable_ph = caltable(gcachep, vec(gmeas_ph))
───────────┬────────────────────────────────────────────────────────────────────
      time │      AA            AP           AZ          JC          LM        ⋯
───────────┼────────────────────────────────────────────────────────────────────
 0.917±0.0 │ 1.0±0.0  -0.814±0.004      missing     missing    -0.5±2.8   2.31 ⋯
 1.217±0.0 │ 1.0±0.0  -0.871±0.003      missing     missing   -0.61±2.8   2.11 ⋯
 1.517±0.0 │ 1.0±0.0  -0.928±0.003      missing     missing   -0.67±2.8   1.91 ⋯
 1.817±0.0 │ 1.0±0.0  -0.978±0.004      missing     missing   -0.75±2.7   1.73 ⋯
 2.117±0.0 │ 1.0±0.0  -1.002±0.004      missing     missing   -0.96±2.7   1.54 ⋯
  2.45±0.0 │ missing       1.0±0.0      missing     missing   2.56±0.33   2.33 ⋯
  2.75±0.0 │ 1.0±0.0  -1.047±0.006      missing     missing    -1.1±2.6   1.12 ⋯
  3.05±0.0 │ 1.0±0.0  -1.133±0.003      missing     missing    -1.1±2.6   0.91 ⋯
  3.35±0.0 │ 1.0±0.0  -1.152±0.004      missing     missing    -1.2±2.6    0.7 ⋯
 3.683±0.0 │ 1.0±0.0  -1.162±0.006    2.36±0.37     missing   -0.88±2.7   0.46 ⋯
 3.983±0.0 │ 1.0±0.0  -1.167±0.003    2.19±0.37     missing   -0.81±2.8   0.24 ⋯
 4.283±0.0 │ 1.0±0.0   -1.16±0.003    2.01±0.37     missing   -0.86±2.8  0.049 ⋯
 4.583±0.0 │ 1.0±0.0  -1.142±0.003    1.81±0.37   -1.7±0.34    0.14±2.9  -0.18 ⋯
 4.917±0.0 │ 1.0±0.0  -1.074±0.006     1.6±0.36  -1.57±0.34   2.65±0.97  -0.42 ⋯
 5.183±0.0 │ 1.0±0.0  -1.079±0.006    1.42±0.36  -1.48±0.34   2.63±0.31  -0.61 ⋯
  5.45±0.0 │ 1.0±0.0   -1.118±0.01    1.27±0.36  -1.37±0.34    2.41±0.3  -0.78 ⋯
     ⋮     │    ⋮          ⋮             ⋮           ⋮           ⋮           ⋮ ⋱
───────────┴────────────────────────────────────────────────────────────────────
                                                    2 columns and 9 rows omitted

Now let's plot the phase curves

plot(ctable_ph, layout=(3,3), size=(600,500))
Example block output

and now the amplitude curves

plot(ctable_am, layout=(3,3), size=(600,500))
Example block output

Finally let's construct some representative image reconstructions.

samples = skymodel.(Ref(post), chain[begin:2:end])
imgs = intensitymap.(samples, fovx, fovy, 128,  128)

mimg = mean(imgs)
simg = std(imgs)
fig = CM.Figure(;resolution=(400, 400))
CM.image(fig[1,1], mimg,
                   axis=(xreversed=true, aspect=1, title="Mean Image"),
                   colormap=:afmhot)
CM.image(fig[1,2], simg./(max.(mimg, 1e-5)),
                   axis=(xreversed=true, aspect=1, title="1/SNR",),
                   colormap=:afmhot)
CM.image(fig[2,1], imgs[1],
                   axis=(xreversed=true, aspect=1,title="Draw 1"),
                   colormap=:afmhot)
CM.image(fig[2,2], imgs[end],
                   axis=(xreversed=true, aspect=1,title="Draw 2"),
                   colormap=:afmhot)
fig

Now let's check the residuals

p = plot();
for s in sample(chain, 10)
    residual!(p, vlbimodel(post, s), dvis)
end
p
Example block output

And viola, you have just finished making a preliminary image and instrument model reconstruction. In reality, you should run the sample step for many more MCMC steps to get a reliable estimate for the reconstructed image and instrument model parameters.

This page was generated using Literate.jl.